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ABSTRACT 


A  model  which  describes  the  thermal-hydraulic  behavior  of  a  packed  particle  bed 
reactor  fuel  element  is  developed  and  compared  to  a  reference  standard  (Tuddenham^c- 
1989).  The  model  represents  a  step  toward  a  thermal-hydraulic  module  for  a  real-time, 
autonomous  reactor  power  controller. 


The  general  configuration  of  the  fuel  element  is  similar  in  construction  to  a  design 
studied  by  Brookhaven  National  Laboratory  and  Sandla^afionat  Laboratory  .^X  bed  of 
small  (diameter  £  500  pm)  fuel  particles  are  packed  between  concentrically  mounted 
retention  cylinders  referred  to  as  frits.  The  element  is  cooled  by  parahydrogen  which 
flows  axially  through  the  inlet  andj>«flet  plenums  and  radially  inward  through  the  fuel 
particle  bed.  r 

The  momentum  integral  approach  used  in  the  MINET  code  (Van  Tuyle,  et  Al, 

1984)  is  applied  to  this  model  to  balance  the  fundamental  mass,  energy  and  momentum 
conservation  relationships.  The  element  is  divided  into  only  three  control  volumes:  the 
inlet  plenum  and  cold  frit  define  the  first  control  volume,  the  fuel  particle  bed  defines  a 
second  control  volume,  and  the  outlet  plenum  and  hot  frit  define  the  third  control  volume. 
The  solid  phase  of  the  particle  bed  is  represented  by  a  single  node. 


This  simple  model  was  validated  against  the  reference  standard  and  compared 
favorably.  As  a  demonstration  of  the  model’s  flexibility,  a  number  of  variations  were 
analyzed.  These  included  variations  in  fuel  element  geometry  and  the  initial  and  final 
values  of  inlet  temperature,  inlet  pressure,  and  outlet  pressure.  As  a  final  demonstration, 
a  cluster  of  nineteen,  1  meter  long  fuel  elements,  arranged  to  form  a  core,  were  analyzed 
for  an  up-power  transient  from  0  MWt  to  approximately  18  MWt. 


The  simple  model  significantly  decreases  the  time  necessary  to  perform  a  single 
analysis.  A  transient  of  10  s  with  a  timestep  of  10  ms,  for  example,  takes  approximately 
45  s  of  computation  on  a  desktop  computer  equipped  with  an  80386  microprocessor. 
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CHAPTER  1 


INTRODUCTION 

1.1  STATEMENT  OF  PROBLEM 

A  packed  particle  bed  reactor  (PBR)  is  a  gas-cooled  reactor  similar  in  nature  to  a 
high  temperature  gas-cooled  reactor  (HTGR).  Unlike  the  HTGR,  however,  the  PBR 
packs  small  fuel  particles  between  inner  and  outer  retention  elements,  designated  as  frits. 
The  PBR  is  appropriate  for  a  number  of  gas-cooled  reactor  applications  and,  in  particular, 
seems  to  be  most  appropriate  for  use  in  space  because  of  its  compactness,  high  outlet 
temperature,  and  wide  range  of  delivered  power. 

The  PBR  is  proposed  for  use  as  a  power  source  for  a  variety  of  systems  in  space. 
The  requirements  of  these  systems  range  from  relatively  low  power  levels,  on  the  order 
of  kilowatts,  to  multi-megawatt  applications  (A-l).  These  systems  include  nuclear  pro¬ 
pulsion  rockets  as  well  as  systems  which  may  require  large,  rapid  power  transients.  The 
proposed  PBR’s  will  initially  be  unmanned  and  will  therefore  require  an  integrated 
autonomous  control  system. 

The  Advanced  Controls  Group  at  MIT  has  worked  with  the  Sandia  National  Labo¬ 
ratory  to  develop  a  control  system  for  a  PBR  which  is  to  operate  in  space.  The  system 
which  is  currently  being  considered  uses  a  control  algorithm  based  upon  a  real-time 
model  which  must  accurately  represent  the  reactor  so  that  proper  control  signals  are  gen¬ 
erated.  As  a  step  toward  the  development  of  the  real-time  model,  a  thermal-hydraulic 
model  was  formulated  by  Tuddenham  (T-l)  to  calculate  the  behavior  of  a  PBR  fuel 
element  during  transient  operations.  This  model  was  developed  to  provide  detailed  infor¬ 
mation  of  fuel  element  behavior  and  to  serve  as  the  standard  for  future  models.  The 
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objective  of  this  investigation  is  to  continue  the  development  of  the  model -based 
algorithm  by  evolving  the  thermal-hydraulic  model  of  a  PBR  fuel  element  and  begin  cou 
pling  the  element  with  other  system  components. 

The  automatic  reactor  controller  will  rely  on  the  thermal-hydraulic  model  for  per¬ 
formance  predictions  based  on  temperature  and  pressure  inputs  from  the  system.  This,  in 
turn,  will  allow  the  controller  to  calculate  projected  fuel  temperatures  and  al’ow  operation 
without  exceeding  safety  limits  within  the  fuel.  Further,  the  model  provides  information 
on  the  temperature  and  pressure  of  the  hydrogen  coolant  which  is  used  to  obtain  the 
amount  of  reactivity  feedback  from  changes  in  coolant  conditions  during  a  transient. 

1.2  APPROACH 

The  major  emphasis  of  this  investigation  is  the  simplification  of  the  thermal- 
hydraulic  reference  model  provided  by  Tuddenham  (T-l).  Although  his  model  produces 
a  very  detailed  analysis,  exceedingly  long  computational  times  are  required  to  reach  a 
solution.  The  simplified  model  allows  for  a  method  of  analysis  which  is  sufficient  in 
detail  to  give  an  accurately  calculation  of  reactivity  feedback  within  the  fuel  and  yet  is 
sufficiently  simple  so  as  not  to  require  excessive  computation. 

Once  the  less  complicated  model  is  developed,  it  must  be  validated  by  comparing 
the  model  results  with  the  reference  standard.  The  reference  standard  provides  calculated 
state  variables  for  some  steady  state  and  transient  conditions. 

In  addition  to  building  a  simplified  model  of  the  fuel  element  behavior,  preliminary 
steps  are  taken  to  couple  the  fuel  element  with  other  fuel  elements  within  the  reactor  as 
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well  as  with  other  core  and  power  plant  components.  Future  developments  can  then  pro¬ 
ceed  with  exact  modeling  of  a  complete  system  once  specifics  become  available  for  each 
separate  component. 

1.3  CONCERNS 

The  major  concerns  involved  with  the  simplification  of  the  reference  standard  are 
those  which  deal  with  the  effects  of  reducing  the  number  of  control  volumes  within  the 
model  and  increasing  the  size  of  the  time  increment.  The  number  of  control  volumes  is 
decreased  from  approximately  fifty  to  only  three.  The  size  of  the  time  step  depends  on 
the  numerical  method  used  and  must  be  sized  to  provide  stability  in  the  formulation.  An 
implicit  method  is  used  in  the  simple  model  and  allows  the  time  step  to  be  increased  from 
.05  ms  to  100  ms  without  loss  of  stability.  Unfortunately,  spatial  resolution  is  decreased 
with  an  increase  in  the  size  of  the  control  volumes  and,  if  spatial  resolution  is  important, 
the  accuracy  of  the  model  will  suffer. 

Specific  control  volume  concerns  include: 

a.  Coolant  Flow:  Can  coolant  flow  through  the  inlet  plenum,  particle  bed, 
outlet  plenum  and  the  retention  elements  be  simplified  into  a  one  dimensional  flow  com¬ 
pared  to  a  two  dimensional  flow  represented  in  the  reference  standard? 

b.  Flow  Distribution:  Can  simplifying  assumptions  be  made  to  use  a  uniform 
flow  distribution  across  the  frits  and  particle  bed? 

c.  Sources  of  Pressure  Loss:  Can  the  sources  of  pressure  loss  in  many  small 
control  volumes  be  combined  in  calculations  for  larger  control  volumes? 

d.  Heat  Transfer:  Can  the  heat  transfer  prediction  of  the  reference  standard 
be  well  represented  by  using  less  control  volumes? 
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Specific  numerical  concerns  include  whether  the  stability  advantages  of  an  implicit 
model  are  sufficient  to  counterbalance  the  longer  computational  times  required. 

The  following  chapters  describe  the  development  and  proposed  application  of  the 
PBR  and  then  describe  the  development  of  a  less  complex  model  to  be  used  to  analyze 
thermal-hydraulic  transient  response  of  a  PBR  fuel  element. 
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CHAPTER  2 


BACKGROUND 


2.1  GENERAL 

The  concept  of  operating  a  nuclear  reactor  system  in  space  is  not  new.  The  require¬ 
ment  to  supply  a  reliable,  sustained  source  of  power  to  satellites  and  deep  space  probes 
was  identified  early  in  the  US  space  program.  Two  programs  were  started  to  investigate 
the  possible  applications  of  nuclear  reactors  in  space.  The  NERVA  (Nuclear  Engine  for 
Rocket  Vehicle  Application)  program  was  started  in  the  early  1960’s  and  was  primarily 
concerned  with  the  development  of  nuclear  technology  for  use  as  a  means  of  propulsion 
in  space.  The  SNAP  (Space  Nuclear  Applications  Program)  program  and  Advanced  Liq¬ 
uid  Metal  Cooled  Reactor  (ALMCR)  program  were  developed  in  parallel  with  NERVA 
and  were  concerned  with  the  development  of  power  supplies  for  satellites  and  deep  space 
probes  (1960  through  1973).  Although  much  of  the  SNAP  program  was  restricted  to  the 
design  and  deployment  of  RTG  (Radio-isotope  Thermal  Generator)  power  sources,  an 
operating  liquid  metal  reactor,  the  SNAP-10A,  was  launched  in  1965  and  subsequently 
tested  at  criticality  in  orbit. 

High  power,  pulsed,  reactor  concepts  are  being  explored  for  orbital  applications  and 
general  boost  applications.  With  the  renewed  emphasis  on  space  exploration,  NTR  (Nu¬ 
clear  Thermal  Rocket)  technology  is  being  reexamined  for  possible  uses  to  minimize 
IMEO  (Initial  Mass  in  Earth  Orbit)  for  long  duration  missions  (L-2).  Although  the 
NERVA  program  has  already  established  an  NTR  technology  base,  other  advanced  con¬ 
cepts  are  being  studied  as  attractive  alternatives  to  this  1960’s  technology. 
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2.2  PARTICLE  BED  REACTORS 


2.2.1  Applications 

Placing  any  payload  into  LEO  (Low  Earth  Orbit)  requires  extensive  planning  and 
support.  This  is  usually  measured  by  the  cost  of  placing  a  unit  weight  into  orbit  (ie,  $/kg) 
and  places  an  emphasis  on  minimizing  IMEO.  Budget  limitations  drive  the  need  for  a 
compact,  lightweight  power  system  which  can  provide  adequate  energy  for  its  applica¬ 
tion.  Conventional  terrestrial  reactor  designs  are  too  bulky  for  use  in  space  so  advanced 
reactor  concepts  are  required.  One  of  these  advanced  concepts  is  the  PBR  which  is 
projected  to  provide  up  to  50%  greater  specific  impulse,  for  a  propulsion  application, 
than  the  current  NERVA  design  (L-2). 

PBR’s  are  being  investigated  by  the  Idaho  National  Engineering  Laboratory  (INEL) 
Multi-Megawatt  Project  Office,  Brookhaven  National  Laboratory  (BNL)  and  Sandia 
National  Laboratory  (SNL)  (H-l,  L-l,  V-l).  Open  cycle  systems  (Figure  2.1)  (M-l)  as 
well  as  closed  cycle  systems  (Figure  2.2)  (G-l)  are  being  examined.  Open  cycle  systems 
are  typical  of  a  PFNTR  (pressure  fed  nuclear  thermal  rocket)  (H-2)  while  closed  cycle 
systems  are  typical  of  pulsed  power  applications  being  studied  by  SNL  within  the  scope 
of  the  PIPE  experiments  being  conducted  in  the  Annular  Core  Research  Reactor  (ACRR) 
(V-l).  PBR  concepts  show  promise  in  that  direct  cooled  particles  within  each  fuel  ele¬ 
ment  provide  a  greater  power  density  and  a  relatively  large  surface  area  for  the  heating  of 
the  working  gas  resulting  in  lighter  and  smaller  systems  (P-1,2). 
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Tank  Pressurization  System 


Figure  2.1:  Typical  Open  Cycle  Nuclear  Thermal  Rocket  System 

(Adapted  from  M-l) 
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Figure  2.2:  Typical  Closed  Cycle  Space  Nuclear  Power  System 
(Adapted  from  G-l) 
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2.2 2  Particle  Bed  Reactor  Design 


The  previous  investigation  of  PBR  thermal-hydraulic  characteristics  was  conducted 
by  Tuddenham  (T-l)  and  addressed  the  complications  that  arise  during  pulsed  operations. 
Pulsed  operations  are  applicable  to  open  as  well  as  closed  cycle  systems.  A  pulse  asso¬ 
ciated  with  the  open  cycle  system  would  be  analogous  to  a  timed  bum  in  a  chemical 
rocket  system. 

The  open  cycle  pulsed  reactor  usually  consists  of  a  coolant  reservoir,  a  coolant 
pump,  a  preheating  stage,  the  reactor  assembly  and  the  exhaust  nozzle.  These  systems 
are  discussed  in  detail  in  references  B-2,  L-2,  P-1  and  P-2.  The  coolant  reservoir  stores 
liquid  hydrogen  which  is  used  as  a  coolant  and  eventually  as  the  exhaust  gas  from  the 
NTR  nozzle.  To  maintain  a  relatively  constant  inlet  pressure  to  the  reactor  assembly  and 
provide  sufficient  mass  flowrate  through  the  core,  a  turbo-pump  is  placed  in  the  system  to 
force  the  coolant  through  the  nozzle  cooling  jacket  and  then  to  the  core.  The  hydrogen 
coolant  is  preheated  becoming  supercritical  as  it  passes  through  the  nozzle  cooling  jacket 
and  innerpass  cooling  channels  within  the  core.  This  provides  efficient  heat  removal 
from  the  nozzle  which  is  subjected  to  high  temperature  exhaust  gases  and  ensures  that  the 
hydrogen  coolant  is  in  the  gaseous  state  prior  to  entering  the  core.  The  hydrogen  then 
passes  through  the  fuel  assemblies,  is  heated,  and  finally  exhausted  into  space  via  the 
nozzle. 
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Figure  23:  Packed  Particle  Bed  Fuel  Element 
(Adapted  from  L-3) 
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As  the  hydrogen  coolant  passes  through  the  fuel  elements,  it  experiences  several 
directional  changes  as  shown  in  Figure  2.3  (L-3).  The  fuel  elements  are  clustered, 
packed  particle  beds  mounted  in  a  moderator  support  assembly  and  surrounded  by  control 
drums  and  reflectors.  A  pressure  vessel,  most  likely  aluminum,  encloses  the  entire 
assembly  (P-2).  The  fuel  elements  may  be  clustered  as  shown  in  Figure  2.4  in  a  number 
necessary  for  the  power  required. 

Each  fuel  element  assembly  consists  of  fuel  particles  packed  between  the  cold  frit, 
the  outer  retention  element,  and  the  hot  frit,  the  inner  retention  element.  These  retention 
elements  are  porous  to  allow  the  coolant  to  flow  through  the  assembly  but  are  fabricated 
to  be  robust  enough  to  retain  the  fuel  particles  between  them  at  elevated  temperatures. 

The  cold  frit  serves  to  distribute  the  coolant  flow  axially  as  well  as  to  retain  the  fuel  par¬ 
ticles.  This  allows  better  cooling  of  more  power  dense  regions  of  the  particle  bed  and 
prevent  preferential  flow  through  specific  regions.  The  cold  frit  is  made  by  sintering 
many  small  (2.5  |im)  stainless  steel  particles  together.  The  hot  frit,  however,  is  made  of  a 
high  temperature  resistant  material,  usually  rhenium,  which  has  evenly  spaced  holes 
drilled  to  attain  a  desired  porosity. 

The  fuel  particles  are  approximately  500  pm  in  diameter  and  consist  of  a  central 
fuel  kernel  surrounded  by  two  pyrographite  layers  and  an  outer  layer  of  zirconium  car¬ 
bide  (Figure  2.5).  Uranium  carbide  is  used  as  the  fuel  and  is  enriched  to  approximately 
93%  for  compactness.  These  particles  are  packed  directly  in  the  annular  region  between 
the  cold  frit  and  the  hot  frit  without  being  imbedded  into  a  graphite  matrix  similar  to  high 
temperature  gas  reactors.  Direct  packing  provides  better  heat  transfer  to  the  coolant  and 
is  expected  to  behave  well  during  rapid  power  transients  associated  with  pulsed  power 
operations  (P-3). 
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2.23  PIPE  Experiment 


Pulsed  Irradiation  of  a  Packed  Bed  Element  (PIPE)  experiments  have  been  con¬ 
ducted  by  Sandia  National  Laboratory  (V-l).  The  PIPE  experiments  attempted  to  evalu¬ 
ate  the  performance  of  a  packed  particle  bed  fuel  element  in  the  areas  of  temperature 
characteristics,  flow  characteristics,  fuel/coolant  interaction,  and  power  output.  A  fuel 
element  test  assembly  (Figures  2.6  and  2.7)  is  then  placed  into  the  annular  core  research 
reactor  (ACRR)  for  evaluation. 

Because  previous  modeling  of  a  PBR  fuel  element  concentrated  on  the  specific 
geometry  of  the  PIPE  experiment,  it  is  important  to  describe  the  experimental  apparatus 
used  for  the  PIPE  experiments.  The  test  assembly  is  a  right  cylinder  4.19  m  long  with  a 
diameter  of  .35  m.  Parahydrogen  is  introduced  to  the  inlet  plenum  of  the  fuel  element  at 
a  pressure  of  2  MPa  and  a  temperature  of  300  K  by  a  series  of  blowers  which  circulate 
the  coolant  through  the  test  assembly.  Flow  enters  the  inlet  plenum  axially  then  flows 
radially  inward  through  the  cold  frit,  the  particle  bed  and  the  hot  frit  The  hydrogen  then 
exits  the  outlet  plenum  into  a  heat  sink  composed  of  stainless  steel  ball,  through  the  flow 
meter  and  back  to  the  blowers  for  another  pass  (V-l). 
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Figure  2.5:  500  Micron  Fuel  Particle 


(Adapted  from  V-l) 
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Figure  2.6:  PIPE  Experiment  Fuel  Element  Test  Assembly 
(Adapted  from  V-l) 
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Figure  2.7:  PIPE  Experiment  Fuel  Element  Test  Subassembly 
(Adapted  from  V-l) 
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Figure  2.8  shows  detailed  dimensions  of  the  fuel  element  tested  by  Sandia  National 
Laboratory.  In  addition  to  the  measurements  shown  in  Figure  2.6,  other  useful  parame¬ 
ters  are: 


•Cold  Frit  Thickness 
•Cold  Frit  Porosity 
•Cold  Frit  Particle  Diameter 
•Cold  Frit  Material 
•Hot  Frit  Thickness 
•Hot  Frit  Porosity 
•Hot  Frit  Material 


1.70  to  2.36  mm 
68.5% 

2.5  pm 

316  Stainless  Steel 
.76  mm 
23.3% 
Rhenium 
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Figure  2.8:  Dimensions  of  PIPE  Experiment  Fuel  Element 
(Adapted  from  T-l) 


2.2.4  Existing  Models 

In  his  thesis  (T-l),  Tuddenham  develops  an  analytical  model  to  calculate  the  behav¬ 
ior  of  the  PIPE  experiment  particle  bed  fuel  element.  His  model,  which  will  serve  as  the 
reference  standard  and  which  will  henceforth  be  referred  to  as  such,  divides  the  fuel 
element  into  50  control  volumes.  The  reference  standard  is  separated  into  10  subdivi¬ 
sions  in  the  radial  direction  and  5  subdivisions  in  the  axial  direction  (Figure  2.9).  Fur¬ 
ther,  it  uses  a  staggered  grid  arrangement  (Figure  2.10)  for  discretization  of  the 
conservation  laws.  Table  2.1  specifies  the  geometry  of  the  reference  standard  and  shows 
how  a  varying  coid  frit  thickness  affects  the  volume  of  the  inlet  plenum  and  the  particle 
bed. 

Although  several  other  models  exist  for  advanced  space  reactors  (B-2,  G-l,  L-2), 
the  reference  standard  addresses  only  the  thermal-hydraulics  of  a  PBR  fuel  element.  It  is 
more  universal  than  one  of  the  models  used  by  Brookhaven  National  Laboratory  (B-l), 
which  assumes  that  power  and  coolant  flow  is  matched  in  each  control  volume.  The  ref¬ 
erence  standard  separates  these  and  provides  detailed  information  pertaining  to  coolant 
temperature,  pressure,  and  density  for  later  use  by  a  neutronic  model.  Also,  the  reference 
standard  attempts  to  focus  on  actual  operating  conditions  by  using  only  instrumented  sys¬ 
tem  parameters,  such  as  inlet  and  outlet  temperatures  and  pressures,  as  sources  of  input. 

Tuddenham  balances  mass,  momentum,  and  energy  in  each  of  his  control  volumes 
at  each  advance  in  t  V  maximum  timestep  of  50  |xs  is  used  by  the  reference  standard. 
Although  very  detailed,  the  calculations  tend  to  take  a  prolonged  length  of  time.  To 
examine  a  6  second  transient,  for  example,  a  single  run  would  take  approximately  48 
hours  of  computational  time  on  an  AT  compatible  desktop  computer  (80286).  On  the 
other  hand,  the  detail  provided  in  the  reference  model  allows  modeling  of  axial  and  radial 
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crossflows  between  the  control  volumes.  Because  the  simple  model  combines  many  of 
the  axial  and  radial  control  volumes  into  a  single  control  volume,  the  particle  bed  for 
instance,  specific  information  concerning  crossflow  is  lost.  If  crossflow  calculations  are 
of  primary  importance,  the  reference  standard  should  be  used. 
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Control  Volume  Definition 

10  Radial  Control  Volumes 
5  Axial  Control  Volumes 
50  Total  Control  Volumes 


Total  Volumes: 
Inlet  Plenum 
Cold  Frit 
Particle  Bed 
Hot  Frit 
Exit  Plenum 


.935 09122  1 
.06269239  1 
.17614176  1 
.00425805  1 
.16905401  1 


Figure  2.9:  Control  Volume  Definition  for  the  Reference  Standard 


Table  2.1:  Reference  Standard  Control  Volume  Geometry 


(measurements  in  mm) 


Axial 

Outer 

Cold  Frit 
Outer 

Cold  Frit 

Outer 

PBED 

Inner  PBED 

Hot  Frit 

Exit 

Plenum 

Position 

Radius 

Radius 

Thickness 

Radius 

Radius 

Thickness 

Radius 

1 

44.58 

29.56 

2.26 

27.30 

15.01 

0.76 

14.250 

2 

44.58 

29.39 

1.93 

27.47 

15.01 

0.76 

14.250 

3 

44.58 

29.33 

1.79 

27.53 

15.01 

0.76 

14.250 

4 

44.58 

29.32 

1.78 

27.54 

15.01 

0.76 

14.250 

5 

44.58 

29.38 

1.91 

27.48 

15.01 

0.76 

14.250 

(Volumes  in  Liters) 


Height 

Volume 

Volume 

Volume 

Volume 

Volume 

Axial 

of  Control 

of  Inlet 

of 

of 

of 

of  Exit 

Position 

Volume 

Plenum 

Cold  Frit 

PBED 

Hot  Frit 

Plenum 

1 

53.00 

981' 

0.021 

0.087 

0.0037 

2 

53.00 

B wM 

0.018 

0.088 

0.0037 

3 

53.00 

0.188 

0.017 

0.089 

0.0037 

0.0338 

4 

53.00 

0.188 

0.017 

0.089 

0.0037 

0.0338 

5 

53.00 

0.187 

0.018 

0.088 

0.0037 

0.0338 

TOTAL: 

265.00 

0.935 

0.092 

0.440 

0.0185 

0.169 

Void 

Fraction: 

1 

0.685 

0.4 

0.23 

1 

Void 

Volume: 

0.935 

0.063 

0.176 

0.0042 

0.169 
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(Adapted  from  T-l) 
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CHAPTER  3 


THE  SIMPLE  MODEL 


3.1  MODEL  BASTS 

A  proposed  particle  bed  reactor  (PBR)  which  is  to  be  placed  in  low  earth  orbit 
(LEO)  will  be  subjected  to  multiple  power  transients.  These  transients  may  take  the  form 
of  start-up  and  shutdown  operations  or  rapid  power  transients  during  pulsed  power  opera¬ 
tions  (B-l,  G-l,  L-l).  Regardless  of  the  form,  any  transient  will  require  precise  control 
of  the  reactor.  As  a  part  of  the  control  system  development,  an  accurate  as  well  as  a 
timely  simulation  must  be  accomplished  via  a  digital  model  of  the  reactor.  When  cou¬ 
pled  with  an  appropriate  neutronic  model,  the  reactor  may  be  safely  controlled  by  analyz¬ 
ing  control  module  input  parameters.  Because  fuel  element  neutronics  are  closely  tied  to 
the  thermal-hydraulic  behavior  of  the  fuel  element,  it  is  important  to  develop  a  good 
thermal-hydraulic  model. 

The  basis  of  this  research  is  to  examine  a  method  of  analysis  for  the  thermal- 
hydraulic  response  of  a  PBR  fuel  element  in  sufficient  detail  as  to  calculate  reactivity 
feedback  within  the  fuel  for  a  subsequent  neutronic  model.  The  analysis  should  be, 
however,  assiduously  simple  so  as  not  to  require  excessive  programming  computations 
experienced  with  the  reference  standard  model.  Additionally,  preliminary  steps  are  to  be 
taken  to  couple  a  single  fuel  element  with  other  fuel  elements  within  the  core. 

The  simple  model  will  have  some  advantages  over  the  reference  standard,  the  most 
important  of  which  is  computational  speed.  Because  the  simple  model  is  faster  than  the 
reference  standard,  many  more  fuel  element  variations  may  be  investigated.  This  permits 
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a  wider  scope  of  transient  analyses  and  allows  a  departure  from  PIPE  geometry.  Further, 
the  output  data  is  presented  quicker,  granting  the  user  a  better  understanding  of  cause- 
and-effect  during  the  fuel  element  design  phase. 

The  following  sections  discuss  the  formulation  of  the  simple  model  and  the  simpli¬ 
fications  made  to  improve  computational  performance  over  the  reference  standard  with 
respect  to  the  amount  of  time  necessary  to  perform  a  single  analysis.  Also  in  this  chapter, 
the  solid  phase  (fuel  particles)  and  the  gas  phase  (coolant)  are  described  in  detail. 

Finally,  a  method  of  solution  is  presented  which  combines  the  gas  and  solid  phases  into 
one  model. 

3.2  SIMPLIFICATIONS 

One  of  the  goals  of  the  simplified  model  is  to  be  able  to  model  the  thermal- 
hydraulic  response  of  a  fuel  element  in  such  a  manner  as  to  develop  a  real-time  analysis. 
A  real-time  model  is  one  which  is  able  to  assess  a  system’s  response  in  a  period  of  time 
less  than  or  equal  to  the  time  it  would  take  for  an  actual  response.  To  do  this,  simplifica¬ 
tions  must  be  made  to  the  reference  standard. 

Merely  loading  the  reference  standard  code  into  a  larger,  faster  computer  (ie,  Cray) 
would  certainly  reduce  the  amount  of  time  necessary  for  a  single  ran.  However,  the  cost 
of  analyzing  the  many  variations  would  preclude  the  necessary  translation.  On  the  other 
hand,  simplification  of  the  mathematical  method  of  the  reference  standard  would  produce 
a  convenient  program  which  could  be  executed  on  a  typical  desktop  computer  in  a  rea¬ 
sonable  time  frame.  Further,  a  version  of  the  model  could  be  considered  as  part  of  a 
spacecraft  power  controller. 
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The  two  major  simplifications  incorporated  into  the  simple  model  are  the  reduction 
of  control  volumes  and  the  increasing  of  the  time  step.  Reducing  the  number  of  control 
volumes  provides  a  proportional  decrease  in  the  amount  of  computation  required  and 
increasing  the  time  interval,  reduces  the  number  of  times  these  calculations  must  be  made 
during  a  specified  simulation  time  interval.  The  overall  effect  is  to  reduce  the  net  number 
of  calculations. 

An  example  of  the  effect  of  the  numerical  simplification  is  readily  seen  when  con¬ 
trasting  the  two  models  on  a  similar  computer  (80386).  It  takes  approximately  8  hours  to 
compute  a  6  second  simulated  transient  using  the  reference  standard  with  a  timestep  of 
.05  ms  but  only  about  2  minutes  for  a  10  second  simulation  transient  using  the  simple 
model  with  a  timestep  of  100  ms.  This  includes  a  reduction  in  the  number  of  control  vol¬ 
umes  for  the  element  from  fifty  in  the  reference  standard  to  only  three  in  the  simple 
model.  A  computational  advantage  of  approximately  240  to  1  is  realized  with  the  simple 
model  which  allows  analysis  of  240  variations  in  the  same  time,  on  a  comparable  desktop 
computer,  that  it  takes  to  do  a  single  variation  using  the  reference  standard. 

The  simple  model  uses  only  three  control  volumes  compared  to  the  fifty  defined  by 
the  reference  standard.  The  three  control  volumes  are  exhibited  in  Figures  3.1,  3.2  and 
3.3.  The  first  control  volume  (Figure  3.1)  includes  the  inlet  plenum  and  the  cold  frit. 

The  second  control  volume  includes  only  the  packed  particle  bed  consisting  of  the  fuel 
particles.  Lastly,  the  third  control  volume  includes  the  hot  frit  and  the  exit  plenum. 
Combining  the  cold  and  hot  frits  with  the  inlet  and  exit  plenums  respectively  allows  a 
more  precise  model  of  the  particle  bed  while  defining  conditions  at  the  inlet  and  outlet  of 
the  fuel. 
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Control  Volume  #1: 
Inlet  Plenum  and  Cold  Frit 


Figure  3.1:  Control  Volume  1 
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Figure  3.2:  Control  Volume  2 
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3.3  SOLID  PHASE 

3.3.1  General 

The  solid  phase  of  this  model  represents  the  fuel  particles  packed  between  the  cold 
frit  and  the  hot  frit.  For  simplicity,  the  particle  bed  is  modeled  as  a  single  control  volume 
and,  to  simplify  further,  the  particle  bed  is  represented  by  a  single,  average  fuel  particle. 
Additionally,  this  average  fuel  particle  is  centered  axially  and  radially  in  the  particle  bed. 
This  approach  is  similar  in  nature  to  the  one  Tuddenham  used  for  each  control  volume  of 
the  reference  standard  and  assumes  that  all  of  the  fuel  particles  behave  in  a  analogous 
manner. 

Fuel  particle  dimensions  and  physical  properties  similar  to  the  PIPE  experiments 
are  used  in  the  modeling  of  the  solid  phase  as  shown  in  Figure  3.4.  Once  the  physical 
properties  of  the  fuel  particle  are  defined,  an  overall  heat  transfer  coefficient  is  developed 
from  a  weighted  average  of  these  properties  and  a  surface  fuel  temperature  may  be  deter¬ 
mined  (M-2). 

Heat  deposition  rate  also  depends  on  the  physical  properties  of  the  fuel  particles 
used  in  the  model.  For  a  power  transient,  the  particle  bed  power  density  is  increased  over 
a  specific  transient  time  to  a  predetermined  level.  The  heat  energy  is  not  immediately 
transferred  to  the  coolant  however.  The  effects  of  heat  deposition,  heat  storage,  and  heat 
removal  must  be  taken  into  account. 
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Figure  3.4:  Fuel  Particle  Geometry 
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3.3.2  Energy  Balance 


The  fuel  particles  are  the  source  of  heat  to  the  hydrogen  coolant  and  can  be  modeled 
as  a  separate  control  volume  from  the  coolant.  Since  ail  of  the  fuel  is  contained  in  a 
single  control  volume,  the  following  equation  for  the  conservation  of  energy  may  be 
employed  for  the  fuel: 


mC  —  =q-hAvVcv(Ts-Tc ) 

V  at  Jcv 

T  =  effective  fuel  particle  temperature 

Ts  =  fuel  particle  surface  temperature 

q  =  heat  source 

h  =  heat  transfer  coefficient 

Av  =  particle  surface  area  per  unit  volume 

Vcv  =  size  of  the  control  volume 

Tc  =  bulk  temperature  of  coolant 


Modeling  of  the  solid  phase  for  this  model  parallels  the  reference  standard.  Refer¬ 
ence  M-2  uses  a  single  node  analysis  for  the  particle  bed.  This  allows  the  use  of  an 
assumed  steady  state  temperature  distribution  through  the  particle  bed.  In  addition, 
material  properties  are  assumed  to  remain  constant  Reference  M-2  simplifies  the  general 
heat  balance  equation  in  3.1  by  developing  weighted  average  value: 


-  dT  -  -  3.2 

maCp  —  =  q”-U(T-Tc) 
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ma  =  particle  mass  per  unit  surface  area  (kg/irf) 

C  p  =  average  particle  specific  heat  (J/kgK) 
q  ”  =  heat  deposition  per  surface  area  (W/m2) 

U  =  effective  over  all  heat  transfer  coefficient  (W/m2K) 


It  must  be  pointed  out  that  T  does  not  represent  a  temperature  at  a  specific  point  in 


the  fuel  particle  but  is  an  average  temperature  for  use  in  equation  3.2.  The  remaining 
values  are  determined  in  accordance  with  the  procedures  of  M-2  which  uses  i  =  1  ...  4  to 
designated  the  shells  of  figure  3.4: 
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Although  the  properties  associated  with  the  fuel  particle  vary  with  temperature,  the 
simple  model,  as  well  as  the  reference  standard,  hold  these  material  properties  constant. 
The  exception  is  the  heat  transfer  coefficient,  h ,  which  is  allowed  to  vary  with  tempera¬ 
ture  and  with  coolant  properties,  h  is  shown  in  the  following  equation  as  a  function  of 
the  Nusselt  number  and  inversely  with  the  sampling  position  in  the  particle  bed.  This 


relationship  will  be  discussed  in  greater  detail  in  the  next  section. 


Koolan,=N“o 


where:  x  =  position  in  the  particle  bed  (in  this  case  x  will  be  half  of  the  particle  bed 
thickness) 

Conductive  heat  transfer  between  adjacent  fuel  particles  and  control  volume  bound¬ 
aries  as  well  as  radiative  heat  transfer,  which  can  occur  at  elevated  temperatures,  are 
neglected  in  the  simple  model.  These  forms  of  heat  transfer  are  addressed  and  discussed 
at  length  in  reference  T-l  which  concludes  that  their  effects  are  minor  compared  to  the 
convective  heat  transfer. 


3.3.3  Discretized  Form 

In  order  to  analyze  the  thermal-hydraulic  response  of  the  PBR  fuel  element,  numer¬ 
ical  approximations  must  be  made  to  the  ordinary  differential  equations  and  the  convec¬ 
tive  equations  used  in  the  model.  The  solid  phase  heat  balance  equation  addressed  in 
equation  3.2  may  be  solved  by  using  a  discretized  form  of  the  equation. 

As  with  the  reference  standard,  each  side  of  equation  3.2  is  multiplied  by  the  con¬ 
trol  volume  size,  V^,  and  the  particle  surface  area  per  unit  volume,  Av,  to  produce  total 
power  terms  (Watts).  This  is  shown  in  the  following  expression: 


3.5 


maCpVcvAv^  =  q"VcvAv  -  UVcvAv(T  -  TG) 


where  Av  can  be  expressed  as  a  function  of  the  void  fraction  and  the  fuel  particle 
diameter,  DP: 

.  _  6(1  -£)  3.6 

To  simplify  this  expression,  incorporate  equations  3.7,  3.8,  and  3.9: 


Let  A  =  UVfyAv 

3.7 

Let  B  =  maC pV (yAv 

3.8 

Let  Q  =q”VcvAv 

3.9 

These  substitutions  in  equation  3.5  result  in: 


dT  -  3.10 

B  —  =  Q-A(T-TC ) 


A  simple  implicit  scheme  is  used  to  place  equation  3.10  into  a  discretized  form. 

This  helps  to  provide  stability  at  larger  timesteps.  In  addition  to  using  T  at  time  n+1 ,  Q 
is  also  advanced  to  n+1.  Qn  +  '  represents  the  driving  element  in  this  relationship.  If  Q 
were  to  remain  unchanged,  the  result  would  be  a  null  transient.  Further,  because  the  A 
(equation  3.7)  is  a  function  of  properties  which  change  with  temperature,  it  is  evaluated  at 
time  n  for  the  discretized  form  of  3.10  which  can  be  expressed  as: 
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5 1 


3.11 


When  this  is  solved  for  Tn+\  a  simple  expression  results: 

r  +  ]  fn+\{Qn  +  '+AnTnG)  3A2 

A* 

1+8 '7 

Equation  3.14  is  incorporated  into  the  simple  model  to  advance  the  average  temper¬ 
ature  of  the  fuel.  The  inteiphase  heat  energy  transferred  from  the  solid  phase  to  the  gas 
phase  can  now  be  calculated  from  the  temperature  difference  between  the  average  fuel 
temperature  and  the  bulk  temperature  of  die  coolant.  The  interphase  heat  transfer,  Q,,  is 
expressed  as: 


Ql"*'  =  VAvvJt,t'-Ti) 


3.13 
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3.4  GAS  PHASE 


3.4.1  General 

Although  less  complicated  in  nature  than  the  reference  standard,  many  of  the  fea¬ 
tures  incorporated  in  modeling  the  gas  phase  in  the  reference  standard  are  included  in  the 
simple  model.  The  gas  phase  is,  however,  more  complex  than  the  solid  phase  in  that  it 
includes  both  hydraulic  and  thermal  characteristics.  Complications  arise,  however,  in 
defining  the  thermal-hydraulic  balances  within  each  control  volume. 

The  general  structure  of  the  gas  phase  equations  balances  energy,  mass,  and 
momentum  for  each  control  volume.  As  enthalpy  and  mass  flowrate  are  solved  for  each 
control  volume,  these  values  are  advanced  in  time  and  forwarded  to  the  adjoining  control 
volume.  This  allows  the  model  to  generate  a  time  response  for  the  fuel  element  during  a 
simulated  transient. 

With  the  exception  of  the  total  number  of  control  volumes,  the  simple  model  is 
based  on  the  same  assumptions  and  significant  features  which  affect  the  element  response 
addressed  in  the  reference  standard.  As  in  the  reference  standard,  the  simple  model  is 
cylindrical  and  includes  axial  and  radial  flows.  Cross  axial  flow  in  the  particle  bed  is  not 
a  variable  for  the  simple  model  since  the  flow  vectors  are  contained  entirely  within  the 
second  control  volume. 

In  the  previous  section,  it  was  mentioned  that  the  heat  transfer  coefficient,  h,  for  the 
hydrogen  coolant  varies  with  temperature  and  velocity.  Reference  E-l  provides  an 
expression  for  the  heat  transfer  coefficient  in  terms  of  the  Nusselt  number,  the  heat  con¬ 
ductivity  of  hydrogen  and  the  sample  position  within  the  particle  bed: 
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where:  Nup  =  Nusselt  number, 


k  =  coolant  heat  conductivity; 
x  =  .5  (particle  bed  thickness) 

Reference  E-l  also  provides  an  expression  for  the  Nusselt  number 

Nup  =  0Me1  Pr33  315 

where:  ovd 

Re=?—£ 


c  _  6<1  ~e> 

r* 


Dp  =  particle  diameter; 

=  effective  particle  diameter  (E-l) 

The  rest  of  this  section  discusses  the  energy,  mass,  and  momentum  balance  tech¬ 
niques  and  equations  used  in  the  simple  model.  Once  the  basic  formulations  are 
identified,  a  discretized  form  of  the  relationships  will  be  stated  and  terms  specific  to  each 
control  volume  will  be  assembled  for  final  analysis. 


3.4.2  Energy  and  Mass  Balance 


The  primary  energy  transfer  mechanism  within  the  particle  bed  is  convection  heat 
transfer  and  it  is  assumed  that  all  of  the  heat  deposited  in  the  fuel  particles  will  be  trans¬ 
ferred  to  the  hydrogen  coolant.  In  order  to  accurately  calculate  the  response  of  the  entire 
fuel  element  to  changes  in  input  heat  energy,  an  energy  balance  and  mass  balance  must 
be  performed  on  each  control  volume.  Since  there  are  fewer  control  volumes  in  the  sim¬ 
ple  model  as  compared  to  the  reference  standard,  the  calculations  involved  will  be  fewer 
as  well. 

The  interphase  heat  transfer  provided  by  equation  3.13  represents  the  heat  energy 
being  convected  to  the  hydrogen  coolant.  This  serves  to  couple  the  solid  phase  with  the 
gas  phase  of  the  model.  The  rate  at  which  heat  is  transferred  from  the  solid  phase  to  the 
gas  phase  depends  on  the  rate  of  increase  of  fuel  heat  deposition  and  coolant  properties. 

For  a  control  volume  with  fixed  volume,  the  time  rate  of  change  of  enthalpy  is  the 
sum  of  the  heat  energy  source,  the  time  rate  of  change  of  control  volume  pressure  and 
volume,  and  the  sum  of  enthalpy  flux  terms.  This  is  similar  to  equation  3.22  in  reference 
T-l: 


(  d(MH) 

{  d, 


=  Q,+v„^-+iwfl, 

at  /  =  i 


3.16 


where:  M  =  mass  within  the  control  volume; 

Qi  =  heat  transfer  from  solid  to  coolant  for  control  volume; 
P  =  average  pressure  within  the  control  volume; 

H  =  specific  enthalpy; 

W  =  mass  flowrate. 
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For  the  simple  model  fewer  control  volumes  limit  the  number  of  enthalpy  flux 


terms  and,  as  a  result,  equation  3.16  may  be  written  as: 


d(MH) 

dt 


3.17 


The  general  mass  balance  equation  may  be  expressed  as: 


—  =  W  -W 
dt 


3.18 


Equations  3.17  and  3.18  may  be  combined  to  produce  a  more  direct  equation  which 
will  allow  us  to  advance  enthalpy.  This  is  accomplished  by  expanding  the  left  hand  side 
of  equation  3.17,  multiplying  both  sides  of  equation  3.18  by  enthalpy,  H,  and  subtracting 
equation  3.18  from  3.17.  Also  assuming  that  Hout  is  representative  of  the  enthalpy  in  the 
control  volume  (a  donor  cell  method)  and  mass  flowrate,  W,  is  the  same  throughout  the 
control  volume,  a  general  relationship  is  produced: 


dHout  dP 

M~dT=Vcv~di+Ql+ W{Hin  ~H°J 


3.19 


In  the  discretized  form  of  equation  3.19,  and  adopting  an  implicit  formulation, 
values  at  the  current  timestep  may  be  used  to  compute  the  enthalpy  for  the  subsequent 
timestep.  An  implicit  numerical  formulation  is  also  used  with  the  enthalpy  to  provide 
numerical  stability  at  large  timesteps.  This  results  in: 
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3.20 


/q 

5 1 


Mn 


n  + 1 


+ 


wn(H r1 


Values  of  {2" +1  and  //"„+1  are  determined  from  the  previous  control  volume  (or  inlet 

boundary)  and  may  be  evaluated  at  n  + 1 .  For  example,  H”n+ 1  for  the  first  control  volume 
represents  the  enthalpy  of  the  coolant  being  supplied  to  the  control  volume.  The  source 
code  allows  a  change  in  the  inlet  temperature  during  a  transient  which  would  subse¬ 
quently  affect  the  enthalpy  of  the  coolant  leaving  the  control  volume. 


3.4.3  Momentum  Balance 


In  addition  to  a  balance  of  mass  and  energy,  a  balance  of  momentum  must  be  dove¬ 
tailed  into  the  response  of  the  fuel  element.  For  the  simple  model,  the  framework  of  the 
MINET  (Momentum  Integral  Network)  code  (V-2)  is  used  to  determine  the  effects  of 
momentum  on  the  flow  of  coolant.  The  MINET  code  develops  equivalent  forces  which 
act  to  resist  the  flow  of  the  coolant  through  each  control  volume.  These  resistive  forces 
can  be  divided  into  the  significant  mechanisms  which  contribute  to  the  head  loss  through 
a  fuel  element.  Unlike  the  MINET  code,  which  accounts  for  spatial  changes  in  mass 
flowrate,  the  simple  model  assumes  that  mass  flowrate  is  the  same  everywhere  (at  any 
time  t). 

The  basic  momentum  integral  equation  (V -2)  used  by  the  simple  model  is 


dW 


Icv-^  =  p;n-po.-UF^i 


3.21 
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where:  lc v  =  L/A  (control  volume  inertia) 


L  =  average  length  of  travel  for  each  control  volume 
A  =  cross  sectional  area  perpendicular  to  flow 
F ^  =  equivalent  resistive  pressure  drop 


This  equation  is  equivalent  to  a  constitutive  relationship  which  relates  a  change  in 
pressure  to  the  flowrate  through  a  control  volume  (C-l): 

3" 


where:  RTotal  =  total  equivalent  resistance; 
AP  =  pressure  drop. 


Equation  3.21  may  now  be  expressed  as  a  function  of  the  total  pressure  drop  across 

the  control  volume  and  the  sum  of  flow  resistances: 

dW  "  ,  3.23 

Icv-x=*Pcv-URJW2 


i  -1 


This  is  a  convenient  expression  because  all  of  the  terms  may  be  determined  from 
existing  relationships.  These  include  relationships  for  an  equivalent  pressure  drop  due  to 
friction,  spatial  acceleration,  expansion  or  contraction,  manifold  mixing,  and  packed 
particle  beds.  Although  the  MINET  code  includes  gravity  terms,  they  are  considered 
negligible.  Expressions  which  differ  from  the  reference  standard  will  be  examined  more 
closely. 
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The  equivalent  pressure  drop  due  to  friction  is  calculated  for  the  inlet  and  the  outlet 


plenums  only  and  takes  the  form  of  equation  3.24: 


where:  Ff  =  equivalent  pressure  drop  due  to  friction; 


Lp  =  plenum  half-length; 
ff  =  0.l38Re~°'151 ; 


D 


“t 


3.24 


3.25 


Since  frictional  pressure  drop  is  directly  related  to  the  distance  that  the  coolant  must 
travel  through  the  control  volume,  an  average  path  length  is  defined  for  the  inlet  and 
outlet  plenum  which  is  half  of  the  overall  element  length.  The  cross-sectional  area  used 
in  this  calculation  is  the  area  in  the  plane  perpendicular  to  the  axial  axis  of  the  fuel 
element.  This  area  is  used  because  the  path  length  in  the  axial  direction  (for  the  inlet  and 
outlet  plenums)  is  so  much  greater  than  the  radial  path  length.  Density  is  the  average  of 
the  density  calculated  at  the  inlet  and  outlet  of  each  plenum. 

Equation  3.25  relates  a  correlation  line  which  was  fit  to  the  log-log  Moody  plot  of 
friction  factors  as  a  function  of  Reynold’s  number,  Re,  and  roughness  in  the  channel 
caust  jy  the  frits  (T-l).  This  improves  on  the  relationships  developed  by  McAdams  and 
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Blasius  to  correlate  friction  factors  for  smooth  pipe  and  is  the  relationship  used  by  the 
reference  standard.  Resistance  due  to  friction  is  determined  in  the  inlet  plenum  and  outlet 
plenum  only  since  the  Ergun  (E-2)  relation  includes  this  effect  for  the  particle  bed. 

Spatial  acceleration  in  a  duct  which  allows  cross  sectional  area  to  change  while 
holding  coolant  density  constant  will  produce  a  change  in  pressure  across  a  control  vol¬ 
ume.  These  conditions  are  experienced  in  the  inlet  and  outlet  plenums  and  is: 

(  1  1  W2  3.26 

Ffl"U2  A-2J2P 

where:  Fa  =  equivalent  pressure  drop  due  to  acceleration. 

The  cross-sectional  areas  in  this  case  are  the  areas  associated  with  the  inlet  and 
outlet  of  the  plenum  control  volumes.  For  example,  for  the  outlet  plenum,  A„  is  the  area 
at  the  exit  of  the  fuel  element  and  Aj  is  the  cross-sectional  area  at  the  inlet  of  the  hot  frit. 
A„  for  in  the  inlet  plenum  is  the  cross-sectional  area  at  the  outlet  of  the  cold  frit  and  A*  is 
the  area  at  the  inlet  to  the  fuel  element.  As  in  the  calculation  for  the  frictional  pressure 
drop,  the  density  used  in  equation  3.26  is  the  average  of  the  density  at  the  inlet  and  outlet 
of  the  control  volume  in  question. 

Equation  3.26  can  not  be  used  for  the  control  volume  containing  the  fuel  particles 
because  the  coolant  experiences  a  change  in  duct  area  and  density.  A  general  form  of  the 
change  in  pressure  caused  by  spatial  acceleration  in  a  duct  which  allows  cross  sectional 
area  and  density  to  change  is  addressed  in  Appendix  A.  This  case  describes  the  changes 
to  the  hydrogen  coolant  as  it  travels  radially  through  the  particle  bed.  The  coolant  will 
experience  a  change  in  cross  sectional  flow  area  as  well  as  a  change  in  density  due  to  the 
transfer  of  heat.  This  takes  the  form  of: 
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3.27 


(FJ 


PBED 


f  A  i  +A2  ' 

2A  jA2  j 


Wv2- 


^  Ai  +A;  ^ 

2A,^2  y 


Wv, 


where:  (Fa)PBED  =  spatial  acceleration  losses  in  the  particle  bed; 


A,  =  cross-sectional  flow  area  at  the  inlet  of  the  particle  bed; 

A2  =  cross-sectional  flow  area  at  the  outlet  of  the  particle  bed; 
v,  =  coolant  velocity  at  the  inlet  of  the  particle  bed; 
v2  =  coolant  velocity  at  the  outlet  of  the  particle  bed. 

Another  source  of  resistance  within  the  fuel  element  is  the  sudden  expansions  and 
contractions  experienced  by  the  coolant  as  it  passes  through  the  cold  frit  and  hot  frit 
These  expansion/contractions  would  occur  when  the  coolant  enters  and  exits  each  frit 
because  of  the  area  change  on  both  sides  of  the  frits.  For  example,  the  coolant  would 
experience  a  first  contraction  as  it  entered  the  cold  frit  and  a  second  contraction  as  it 
leaves  the  cold  frit  and  enters  the  particle  bed.  The  coolant  would  experience  a  first 
expansion  as  it  leaves  the  particle  bed  and  passes  into  the  hot  frit  and  a  second  expansion 
as  it  enters  the  outlet  plenum.  The  general  form  of  this  equivalent  pressure  drop  is 
adopted  from  T-2  and  is  given  as: 


Fk  =  Kk 


W 2 


2pA 


small 


3.28 


where:  Fk  =  equivalent  loss  due  to  a  sudden  expansion  or  contraction; 


Kk  =  K.  for  expansions; 
*,  =  d-B)2; 

Kk  =  Kc  for  contractions; 
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g  smaller  area 
larger  area  ’ 

A-snaii  =  the  smaller  of  the  two  areas  involved. 

The  value  of  density  used  in  equation  3.28  is  the  density  at  the  point  of  expansion  or 
contraction.  In  the  case  of  the  cold  frit,  the  density  calculated  at  the  entrance  of  the 
particle  bed  is  used  to  evaluate  the  equivalent  pressure  drop  due  to  sudden  contractions. 
The  density  calculated  at  the  exit  of  the  particle  bed,  on  the  other  hand,  is  used  for  the 
expansions  observed  at  the  hot  frit. 

The  final  relationship  to  be  used  to  describe  the  hydraulics  of  the  fuel  element  is  for 
the  mixing  effects  experienced  in  the  inlet  and  outlet  plenums.  This  is  known  as  man¬ 
ifold  mixing.  The  general  relationship  was  derived  by  Bajura  (B-4,  5, 6)  and  improved 
by  Datta  and  Majumdar  (D-l,  M-3),  and  is  used  in  the  reference  standard.  When  put  into 
a  form  compatible  with  equation  3.22,  the  manifold  effect  is  expressed  as: 


Fm  —  -  radial  +  , 


M  —  axial 


W2 

PALo/ 


3.29 


M  — radial 


W2 

P^radial 


where:  Fm  =  equivalent  pressure  drop  due  to  the  manifold  effect 


0  =  .95  for  the  inlet  plenum  and  1 . 1  for  the  exit  plenum 


The  value  of  density  used  in  this  case  is,  again,  the  average  of  the  inlet  and  outlet 
densities  for  the  control  volume  being  analyzed.  In  addition,  it  must  be  noted  that 
although  the  reference  standard  uses  a  value  of  2.66  for  exit  plenum  0  in  the  calculations. 
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T-l  states  that  the  values  of  0  may  be  adjusted  to  take  into  account  differences  in  mixing. 
Because  of  the  difference  in  how  the  simple  model  treats  crossflow  within  a  control 
volume,  0  for  the  outlet  plenum  is  reduced  from  2.66  to  1.1.  This  allows  a  better 
approximation  of  the  results  obtained  by  the  reference  standard.  Although  this  is  minor 
contributor  to  the  total  flow  resistance,  it  serves  to  illustrate  how  the  model  may  be  fine 
tuned  and  results  in  a  better  agreement  for  pressure  drop-flow  relations. 

To  complete  the  set  of  expression  used  in  the  simple  model,  the  Ergun  relation  is 
used  to  determine  the  equivalent  pressure  drop  across  the  cold  frit  and  the  particle  bed. 
Since  the  cold  frit  is  made  of  small  particles,  modeling  it  as  a  particle  bed  is  a  good 
approximation.  The  Ergun  relation  treats  the  flow  through  a  packed  particle  bed  as  flow 
through  tubes  with  irregular  cross  sections  vice  flow  around  many  objects  and  is  the  gen¬ 
erally  accepted  approach  to  modeling  particle  beds  (B-6,  E-2).  The  general  form  of  the 
Ergun  relation  is: 


j  150p.vo(l  -e)2|  f  1.75pv0 1  v0 1  (l  -e) 

M  Die  PM  D„e  j 

where:  Fpbtd  =  equivalent  pressure  drop  due  particle  bed  effects; 

Lb  =  thickness  of  the  particle  bed; 

Dp  =  particle  diameter; 
v0  =  superficial  velocity; 

(J.  =  viscosity; 

e  =  void  fraction  or  porosity; 
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Using  the  conservation  of  mass  relation,  this  expression  may  be  put  into  a  form 


which  is  compatible  with  equation  3.27  and  will  appear  as: 


FptxJ  ~ 


L 

f  I50p(l  -e)2> 

+u 

(  1.75(1  —  E)  Y| 

1  ‘ 
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b 

l  DpeP At  Jj 

3.31 


where: 


W 

PV0 


Once  more,  the  density  used  in  this  relationship  an  average  of  the  density  at  the 
particle  bed  inlet  and  the  density  at  the  exit  of  the  particle  bed. 

Having  defined  all  of  the  significant  terms  used  to  describe  the  pressure  loss 
through  the  fuel  element,  the  next  step  is  to  discretize  equation  3.21.  An  implicit  scheme 
may  be  used  if  W2  is  separated  into  a  combination  of  the  mass  flowrate  evaluated  at  n  and 
n+1  to  give:WnW"  +  '.  The  discretized  form  of  equation  3.21  is  then: 


(W 


n  +  1 


Wn)  1  .  .  „ 

Sr  ~/cv(  cv  roM' 


WnWn  +  ') 
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Solving  for  Wn  + ’: 


Wm  +  t 


WnIcy  +  btAPZv 
Icy  +  StR^W- 


3.33 
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This  section  has  described  the  process  by  which  the  simple  model  is  able  to  advance 
enthalpy  and  mass  flowrate.  The  following  section  will  describe  how  this  is  coupled  with 
the  advance  of  fuel  temperature  and  interphase  heat  transfer  to  produce  a  full  portrayal  of 
the  thermal-hydraulic  response  of  a  typical  fuel  element. 
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3.5  METHOD  OF  SOLUTION 


A  simple  flow  chart  of  the  transient  solution  is  shown  in  Figure  3.5.  As  can  be  seen 
in  the  schematic,  initial  conditions  are  established  by  the  user  and  an  initial,  steady  state 
condition  is  calculated.  This  provides  initial  state  variables  to  the  transient  portion  of  the 
source  code.  Next,  enthalpy  and  mass  flowrate  is  advanced  for  each  control  volume.  In 
addition,  in  control  volume  2,  calculations  include  the  advance  of  fuel  temperature.  This 
in  turn  provides  the  calculated  value  of  the  interphase  heat  transfer. 

Of  note,  the  control  variables  (heat  deposition,  inlet  pressure,  outlet  pressure,  and 
inlet  temperature)  are  also  advanced  in  accordance  with  the  transient  parameters  selected 
at  the  initiation  of  the  simulation. 

Coolant  properties,  such  as  density,  viscosity  and  heat  transfer  coefficient,  are  recal¬ 
culated  for  each  control  volume  based  on  new  temperatures  and  pressures  corresponding 
to  the  value  of  enthalpy  and  mass  flowrate.  Finally,  fuel  element  information  is  sent  to 
the  output  file  and  time  is  advanced  by  the  chosen  timestep  for  the  next  set  of  calcula¬ 
tions. 
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INPUTS 


INITIAL  VALUES 


Heat  Deposition  (q’”) 
Inlet  Pressure  (Pin) 
Outlet  Pressure  (Pout) 
Inlet  Temperature  (Tin) 


Calculations  to  Establish 
Initial  Conditions 


ADVANCE  CONTROL  VARIABLES  - - 

(Pin,  Pout,  Tin,  q’”) 

Control  Volume  1 

ADVANCE  H 

ADVANCE  W  & 

K 

V 

Control  Volume  2  g 

ADVANCE  Tf  P 

ADVANCE  H  * 

ADVANCE  W  S 

o 

«-> 

Control  Volume  3  » 

ADVANCE  H  S 

ADVANCE  W  i 

EVALUATE  PROPERTIES  OF 
COOLANT  AT  NEW  T,  P,  W  - 

I 

OUTPUT  TO  DATA  FILE 


Figure  3.5:  Flow  Chart  for  Transient  Solution 
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CHAPTER  4 


VALIDATION  AND  APPLICATIONS 

4.1  GENERAL 

Although  originally  intended  as  an  interim  step  toward  the  development  of  a  real¬ 
time  controller  for  a  particle  bed  reactor,  the  simple  model  also  lends  itself  serendipi- 
tously  as  a  design  tool.  Because  the  time  necessary  for  a  single  analysis  is  significantly 
shorter  than  that  required  by  the  reference  standard,  many  more  variations  may  be 
examined  and  compared.  Prior  to  using  the  source  code  to  examine  these  different  varia¬ 
tions,  the  source  code  must  first  be  verified  as  a  valid  model.  This  may  be  done  by  com¬ 
paring  the  results  from  the  simple  model  with  the  reference  standard. 

An  additional  feature  of  the  simple  model  is  that  it  is  capable  of  analyzing  transients 
in  which  control  variables  change  (inlet  pressure,  outlet  pressure,  inlet  temperature,  and 
heat  deposition  rate).  That  is,  inlet  pressure,  for  example,  may  be  permitted  to  increase  or 
decrease  during  the  transient. 

4.2  VALIDATION 

Validation  of  the  simple  model  is  accomplished  by  comparing  results  of  steady  state 
and  transient  calculations  provided  by  each  model.  Also,  as  a  verification  that  the  numer¬ 
ical  scheme  is  valid,  null  transients  are  analyzed  at  three  different  power  density  levels. 
Overall,  the  simple  model  agrees  very  closely  with  the  reference  standard  but  lags 
slightly  during  a  baseline  1  second  transient  from  0  GW/m3  to  2  GW/m3.  Null  transients 
are  performed  at  0, 1,  and  2  GW/m3. 
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A  more  in  depth  discussion  of  the  validation  procedure  is  attached  as  Appendix  C 
and  includes  graphical  comparisons  of  the  simple  model  and  the  reference  standard. 

After  completing  a  successful  validation,  the  simple  model  may  be  used  to  examine 
variations  other  than  the  baseline  (PIPE)  configuration  used  by  the  reference  standard  and 
the  PIPE  experiment.  Examples  of  which  are  supplied  in  the  following  sections. 

4.3  DESIGN  VARIATIONS  IN  FUEL  ELEMENT  GEOMETRY 

The  source  code  for  analyzing  the  response  of  the  packed  particle  bed  fuel  element 
allows  changes  in  the  dimensions  of  the  fuel  element.  Geometry  variables  which  may  be 
changed  include  outer  element  radius,  outer  cold  frit  radius,  cold  frit  thickness,  particle 
bed  thickness,  hot  frit  thickness,  and  length.  Other  parameters  which  may  be  changed 
which  influence  the  hydrodynamic  characteristic  of  the  fuel  element  include  particle 
diameter,  frit  porosity,  manifold  mixing  coefficients,  and  fuel  void  fraction. 

As  a  demonstration,  mass  flowrates  for  a  range  of  fuel  element  power  densities 
were  determined  for  a  number  of  different  fuel  element  lengths  and  plotted  against  the 
thermal  power  that  they  would  deliver.  Figure  4.1  shows  the  results  of  these  calculations 
for  lengths  of  .265  m  (baseline),  .500  m,  .750  m,  and  1.000  m.  In  each  case,  inlet  pres¬ 
sure  was  maintained  at  2000  kPa;  outlet  pressure  was  maintained  at  1915  kPa;  and  inlet 
temperature  was  maintained  at  300  K.  Also,  each  curve  has  eleven  points  of  data  which 
represents  calculations  based  on  different  power  densities  ranging  from  0  GW/m3  (ther¬ 
mal  power  equal  zero)  to  1  GW/m3  (the  end  point  of  each  curve)  at  .1  GW/m3  intervals. 

For  an  open-cycle  nuclear  thermal  rocket  application,  Figure  4.1  could  aide  the 
designer  chose  a  range  of  power  operations  given  specific  mass  flowrate  restrictions.  Or, 
on  the  other  hand,  the  designer  could  calculate  the  amount  of  hydrogen  fuel/coolant 
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required  given  power  level  and  fuel  element  size  restraints.  Figure  4.1  could  be  paired 
with  a  plot  of  exit  temperature  verses  thermal  power  to  ensure  component  temperature 
limits  are  not  exceeded. 


Mass  Flowrate  vs  Thermal  Power 
For  Fuel  Elements  of  Various  Lengths 


Mass  Flowrate  (g/s) 


Thermal  Power  (MW) 


Fuel  Element  Length 

L=.265  m 

- * -  L=.500  m 

X  - 

L=.750m 

— 0 —  L=1.00m 

Figure  4.1:  Mass  Flowrate  vs  Power  for  Various  Element  Lengths 
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4.4  VARIATIONS  IN  THE  MANNER  OF  TRANSIENT  CONTROL 

The  code  has  the  provisions  to  fix  the  initial  and  final  values  of  inlet  and  outlet 
pressure  and  inlet  temperature  at  a  single  value.  They  may  also  be  changed  linearly  dur¬ 
ing  a  portion  of  a  transient  of  specified  duration.  The  time  during  which  one  of  these 
variables  may  change  will  typically  start  when  the  power  transient  begins  but  is  not 
coupled  to  the  same  transient  duration.  For  example,  power  density  may  be  increased 
from  0  GW/m3  to  1  GW/m3  in  1  second  and  inlet  pressure  decreased  90  kPa  over  a  period 
of  5  seconds.  To  conserve  fuel,  a  nuclear  thermal  rocket  may  be  operated  in  such  a  fash¬ 
ion,  that  is,  the  reactor  is  started  up  prior  to  admitting  the  hydrogen  coolant  through  the 
fuel  element.  Figures  4.2  and  4.3  illustrate  this  operating  mode  for  a  1  m  long  fuel 
element  and  show  how  mass  flowrate  and  delivered  thermal  power  (the  product  of  mass 
flowrate  and  enthalpy  change)  respond  to  such  a  transient 

Of  note.  Figure  4.2  indicates  that  thermal  power  will  overshoot  the  required  neutron 
power  due  to  the  extended  decrease  in  the  outlet  pressure.  While  an  increase  in  neutron 
power  tends  to  increase  the  outlet  temperature,  a  decrease  in  outlet  pressure  tends  to 
increase  mass  flowrate.  Thermal  power  continues  to  increase  until  outlet  pressure 
reaches  the  required  operating  value  of  1900  kPa. 
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Thermal  Power  Response  to  Combined 
Is  Power/ 1 5s  Outlet  Pressure  Transient 
(L=lm;  P-Inlet=2000  kPa;  T-Inlet=300  K) 


Power  (MW) 


Pressure  (kPa) 


Nuclear  Heat  Gen 

Thermal  Power 

Outlet  Pressure 

Transient  Starts  at  t=ls 

Figure  4.2:  Thermal  Power  Response  to  Combined  Transients 


4.5  CORE  REPRESENTATION  BY  MULTIPLE  FUEL  ELEMENTS 


As  a  further  demonstration  of  how  the  simple  model  may  be  used  as  a  design  tool,  a 
cluster  of  19  fuel  elements  are  arranged  into  three  rings  around  a  center  fuel  element  as 
presented  in  Figure  4.4.  The  center  element  is  designated  as  O  and  is  surrounded  by  the 
A-ring,  B-ring,  and  C-ring.  Using  a  zero  order  Bessel  function  of  the  first  kind  to 
approximate  the  relative  power  densities  of  each  ring,  the  simple  model  was  used  to  cal¬ 
culate  the  time  behavior  of  each  ring  of  the  reactor  for  mass  flowrate,  thermal  power,  and 
exit  temperature.  Length  was  set  to  1  m,  inlet  pressure  to  2000  kPa,  outlet  pressure  to 
1900  kPa,  and  inlet  temperature  to  300  K. 

Figure  4.5  shows  the  mass  flowrate  response  to  a  4  second  power  transient  for  a 
single  fuel  element  in  each  ring.  The  heat  deposition  rate  is  increased  from  zero  to  1 
GW/m3  in  element  O  with  the  other  rings  being  increased  proportionally.  As  expected, 
the  fuel  elements  with  the  lower  power  densities  develop  less  flow  resistance  and,  as  a 
result,  mass  flowrate  is  greater  through  these  elements  at  the  final  power  level.  Mass 
flowrates  specific  to  each  element  may  be  combined  to  describe  the  overall  mass  flowrate 
response  to  the  transient  as  shown  in  Figure  4.6  such  that: 
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wT=iwi 

1-1 


4.1 


Figure  4.7  depicts  total  thermal  power  response  to  the  4  second  increase  in  neutron 
power.  The  total  thermal  power  is  found  by  summing  the  product  of  mass  flowrate  and 
change  in  enthalpy  for  each  element: 
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Finally,  Figure  4.8  describes  exit  temperatures  for  each  ring  of  elements.  The 
average  exit  temperature  represents  the  temperature  of  the  coolant  which  is  delivered  to 
downstream  power  plant  components  and  is  determined  by  weighting  each  exit  tempera¬ 
ture  with  the  corresponding  mass  flowrate.  Note  that  the  center  ring  will  be  operated  well 
above  the  average  outlet  temperature  which  indicates  that  a  limit  associated  with  fuel 
element  material  may  be  reached  before  any  other  component  (ie,  turbine  blades). 

As  a  design  consideration  factor,  outlet  temperatures  may  be  made  more  uniform  by 
changing  the  orificing  of  the  inlet  plenum  or  the  thickness  of  the  cold  frit  for  the  A,  B, 
and  C-rings.  This  would  serve  to  decrease  thermal  gradients  across  regions  of  the  core  as 
well  as  improving  allowable  thermal  power. 
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Mass  Flowrate  Response  tc  4s 
Up-Power  Transient  for  19  Element 
Particle  Bed  Reactor 


Total  Mass  Flowrate  Response  to  4s 
Up- Power  Transient  for  19  Element 
Particle  Bed  Reactor 
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Thermal  Output  Power  Response  to  4s 
Up-Power  Transient  for  19  Element 
Particle  Bed  Reactor 


Power  (MW) 


Thermal  Power  Neutron  Power 

Transient  Starts  at  t=ls 

Figure  4.7:  Total  Power  Response  for  19  Element  PBR 


Outlet  Temperature  Response  to  4s 
Up-Power  Transient  for  19  Element 
Particle  Bed  Reactor 


4.6  REAL-TIME  CONTROLLER 


The  simple  model  takes  a  step  toward  development  of  a  real-time  model  for  incor¬ 
poration  in  control  of  a  space-based  reactor  by  significantly  decreasing  the  time  necessary 
for  computation.  Although  some  information  is  lost  when  compared  to  the  data  delivered 
by  the  reference  standard,  the  overall  framework  which  describes  fuel  element  perform¬ 
ance  is  similar. 

As  a  demonstration  of  real-time  control,  the  simple  model  was  modified  in  a  man¬ 
ner  which  suspended  the  process  of  writing  data  to  a  data  file.  The  calculation  time  for  a 
transient  was  less  than  the  transient  simulated  (for  a  single  element).  A  10  second 
transient,  for  example,  required  approximately  3  seconds  of  computation  time  for  a  time- 
step  of  100  ms  and  approximately  7.5  seconds  for  a  timestep  of  10  ms.  The  calculations 
were  performed  on  a  desktop  computer  equipped  with  a  30386  microprocessor  (16  MHz). 

Although  there  are  reasons  that  these  results  are  inappropriate  for  a  final  control 
application,  they  do  indicate  the  plausibility  of  using  a  version  of  this  thermal-hydraulic 
model.  The  reasons  for  the  continued  adjustment  of  these  results  are: 

•  Improvements  in  computer  speed  are  expected; 

•  An  optimization  of  calculation  techniques  has  not  been  performed;  and 

•  Multiple  elements  and  other  computations  must  be  done  in  parallel. 
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CHAPTER  5 


CONCLUSIONS  AND  RECOMMENDATIONS 

5.1  CONCLUSIONS 

The  previous  sections  describe  an  investigation  to  construct  a  simple  model  which 
could  effectively  analyze  the  thermal-hydraulic  response  of  a  packed  particle  bed  reactor 
fuel  element  during  steady  state  and  transient  conditions.  The  fuel  element  was  first 
divided  into  a  solid  phase  and  a  gas  phase  which  are  modeled  separately  then  coupled  for 
the  final  solution. 

The  solid  phase  is  modeled  in  a  single  control  volume  and  a  lumped  parameter 
approach  is  used.  This  approach  is  essentially  the  same  as  the  approach  used  in  the  refer¬ 
ence  standard  and  generates  a  fuel  temperature  which  is  representative  of  an  average  fuel 
temperature  instead  of  a  local  fuel  centerline  temperature  for  instance.  Interphase  heat 
transfer  is  calculated  from  the  product  of  an  overall  heat  transfer  coefficient  and  the  dif¬ 
ference  between  the  average  fuel  temperature  and  the  bulk  temperature  of  the  coolant. 
This  interphase  heat  transfer  rate  provides  the  heat  input  to  the  coolant  as  it  passes 
through  the  fuel  particles. 

The  gas  phase  is  modeled  in  three  control  volumes:  the  inlet  plenum  combined  with 
the  cold  frit;  the  packed  particle  bed;  and  the  outlet  plenum  combined  with  the  hot  frit. 
This  particular  arrangement  allows  the  fuel  region  to  be  modeled  separately  from  other 
components  of  the  fuel  element. 

Conservation  equations  for  the  conservation  of  energy,  mass,  and  momentum  (in 
the  form  of  the  momentum  integral  equation)  are  used  to  ensure  a  balance  when  analyz- 
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ing  the  fuel  element  during  a  transient.  Steady  state  values  and  transient  values  (for  a  1  s 
baseline  transient  from  0  to  2  GW/m3)  are  analyzed  and  then  compared  to  the  reference 
standard. 

The  calculations  resulting  from  the  simple  model  compare  favorably  with  the  refer¬ 
ence  standard  and,  because  of  the  simplifications,  are  able  to  analyze  the  same  transient 
significantly  faster  as  well  as  compute  steady  state  values  directly.  After  eliminating 
writing  to  a  data  file,  the  simple  model  is  able  to  achieve  a  faster-than-real-time  thermal- 
hydraulic  analysis  of  a  transient  for  timesteps  larger  than  10  ms  and  an  80386  micropro¬ 
cessor. 

The  simple  model  may  also  be  used  as  a  design  tool.  It  may  be  used  on  a  desktop 
computer  and  many  different  calculations  may  be  accomplished  during  a  reasonable 
timeframe.  Further,  the  simple  model  has  the  capability  to  accept  variations  in  geometry, 
fuel  composition,  and  initial/final  conditions. 

5.2  RECOMMENDATIONS  FOR  FURTHER  STUDY 

Development  of  the  simple  model  is  another  step  in  the  creation  of  a  real-time, 
autonomous  controller  for  a  packed  particle  bed  reactor.  Although  this  model  represents 
the  thermal-fly'*  aulic  portion  of  the  controller,  it  sets  the  stage  for  further  investigation. 

Some  areas  which  would  merit  further  investigation  include: 

•  Calculations  obtained  by  the  reference  model  and  the  simple  model  should  be 
compared  to  experimental  results.  At  the  time  this  investigation  concluded,  such  exper¬ 
imental  results  were  not  available. 

•  A  neutonics  module  needs  to  be  developed  which  uses  the  data  generated  by  the 
thermal-hydraulics  module.  A  real-time  controller  may  then  be  constructed  and  tested. 
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•  A  model  which  includes  additional  power  generating  system  components  should 
be  developed.  Components,  such  as  a  high  temperature  turbine,  turbo-pump,  and  asso¬ 
ciated  valves  and  piping  could  be  coupled  to  the  representation  of  the  core  provided  by 
the  simple  model. 

•  Study  a  means  to  adjust  the  simple  model  to  correct  the  observed  sluggishness 
(Section  C.3) 

•  Continue  to  compare  the  results  obtained  by  the  simple  model  with  future  models 
which  may  construed  as  a  reference  standard. 
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APPENDIX  A: 


DERIVATION  OF  THE  GENERAL  CASE  FOR  SPATIAL  ACCEL¬ 
ERATION  IN  A  DUCT  WITH 

NON-CONSTANT  AREA  AND  NON-CONSTANT  DENSITY 


A.l  BACKGROUND 

Spatial  acceleration  term  in  a  momentum  balance  equation  are  usually  presented  as 
the  pressure  change  in  a  duct  caused  by  a  change  in  density  or  by  a  change  in  cross  sec¬ 
tional  area.  In  the  former  case,  the  pressure  change  in  a  duct  caused  by  a  change  in  den¬ 
sity  while  maintaining  cross  sectional  area  constant  is  expressed  as: 

P>-f2  =  PA!-p,vf  A. 

In  the  latter  case,  a  change  in  pressure  associated  with  spatial  acceleration  in  a  duct 
with  constant  density  but  with  a  change  in  duct  area,  is  expressed  as: 


Pl~P 2  = 


_L_-L 

,  A2  A2 
\A2 


W^_ 

2p 


A.2 


where  W  =  Mass  Flowrate 

In  the  case  of  the  packed  particle  bed  reactor  fuel  element,  the  working  gas  is 
subjected  to  a  combination  of  cross  sectional  area  change  and  density  change  as  it  travels 
in  the  radial  direction.  Although  there  are  integration  methods  available  which  could 
accurately  describe  this  fluid  behavior,  it  would  necessitate  dividing  the  particle  bed  into 
many  control  volumes.  For  the  simple  model,  it  is  desirable  to  maintain  the  particle  bed 
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within  a  single  control  volume.  Because  neither  equations  A.l  nor  A.2  accurately 
describe  the  thermal-hydraulics  associated  with  variable  density  gas  through  the  particle 
bed,  a  general  case  derivation  must  be  formulated  for  the  simplified  model. 

A.2  THE  GENERAL  CASE 

Figure  A.  1  shows  the  general  case  in  which  the  cross  sectional  area  in  the  duct  and 
the  density  of  the  gas  is  allowed  to  change  as  the  gas  travels  through  the  duct.  The  objec¬ 
tive  of  this  derivation  is  to  formulate  an  expression  which  will  satiety  the  limiting  cases 
described  in  equations  A.l  and  A.2  above.  To  do  this  we  must  first  define  an  average 
pressure  P  as  a  combination  of  P}  and  P2: 

P=aPI  +  (l-a)P2  A.3 

Where  a  and  (1  -  a)  represent  fractional  constants  between  0  and  1. 
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To  begin  the  derivation,  a  momentum  balance  analysis  must  be  performed  on  the 
control  volume  shown  in  Figure  A.l.  P  becomes  important  at  this  point  because  it  is 
allowed  to  act  at  some  point  within  the  duct  on  an  area  which  is  defined  as  the  difference 
between  A,  and  A2  when  projected  onto  the  vertical  plane.  This  then  results  in  the 
momentum  balance  equation: 

P  |Aj  —  P  jAj  +  P  (A2  —  Aj)  =  P2V2A2  —  PjVj  A]  A.4 

Substituting  the  definition  of  P  into  equation  A.4  and  collecting  terms  gives: 


(Pi  ~p2)  {A,+a(A2-A,)}  =  p2v22A2  -  p,v,2A, 


A.5 


By  inspection,  equation  A.5  reduces  to  equation  A.l  when  A,  and  A2  are  set  equal 
to  one  another.  When  the  duct  inlet  area  and  outlet  areas  remain  the  same,  the  term 
a(A2-Aj)  is  eliminated  resulting  in  the  limiting  case  shown  in  equation  A.l.  For  the 
other  case  in  which  the  solution  is  known,  the  fractional  constant  a  is  solved  while 
density  remains  constant. 

For  the  case  in  which  density  is  held  constant,  the  Bemolli  relationship  is  used  to 
find  an  expression  for  the  pressure  change  on  the  left  hand  side  of  equation  A.5: 


Pi 


A.6 


and,  because  mass  is  conserved  in  a  steady  state  condition,  constant  cis  used  to  demote  a 
ratio  of  areas  and  hence  a  ratio  of  velocities: 


78 


c=^ 

^2 


Substituting  the  expression  for  P,  —P2  and  letting  v2  =  cv,  results  in  the  following 


expression: 


-(c2-l){A1+fl(A2-A1)}=A,(c-l) 

The  fractional  constant  a  can  now  be  solved  by  dividing  both  sides  by  A2(c  -  1): 


(c  +  l){c  +a(l  -c)}  =  2c 
c(c  +  l)  +  <z(c  +  l)(l-c)  =  2c 
a(c  +  l)(c  - 1)  =  c(c  - 1) 


and  upon  completion  of  the  algebra,  an  expression  is  derived  for  the  fractional  constant 


+  A2 


1  -a  =■ 


A  i  +A2 
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Equations  A.  10  and  A.  11  are  then  substituted  into  equation  A.3  for  an  expression 


for  P: 


P  = 


A, 


A,  +A 


■2J 


At  +A‘ 


2) 


A.12 


Finally,  substituting  this  new  expression  for  P  into  equation  A.4  and  solving  for  the 
pressure  change  provides  an  expression  for  the  general  case: 


Pi 


f 

V 


A,  +  A2 
2A,A2  y 


Wv2 


Aj  +A2 
^  2AjA2 


A.13 
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APPENDIX  B: 


SUMMARY  OF  EQUATION  FOR  HYDRAULIC  ANALYSIS 
OF  A  PBR  FUEL  ELEMENT 


I.  Constitutive  Relationship; 

A Pt  =  RtW2 

APt  =  Total  pressure  drop  across  element 

Rt  =  A  term  representing  the  total  resistance  to  fluid  flow 
W  =  Mass  flowrate  through  the  element 


2.  Momentum  Integral  Equation: 

MINET  CODE  (Reference  V-2) 


r 

m  dt 


=  P, 


-P+F, 


Press  Equivalent 


^PressEqui  valent  Z(Ft+Ff+Fa+Fv+Fp+Fk) 

^ ' PressEqui  valent  =  +Ff+Fa  +Py  +Fp  +RkW2 

RT  =  'L(RI;+Rf  +  Ra+Rv+Rp+Rk) 
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3.  Pressure  Relationships  for  Inlet  and  Outlet  Plenums: 


f  i  \ 


Ff  =  ~ff 


\D«iJ 


W 2 


2p A7 


ff=0.mRe 


-0.151 


L  =  The  half-length  of  the  plenum 


l_  J_]W2 
A?  J2p 

w2 


2pA 


small 


Kk  =  KJor  expansions 
Kk  =  KJor  contractions 


Kc=^(l-Bf 

_  smaller  area 
largerarea 


Fm  —  ~®^Fm -radial  +  -axial) 


M  —  axial 


P^lxial 


and 


M  -  radial 


W2 

P^radial 


0  =  .95  for  Inlet  Plenum;  1.10  for  Outlet  Plenum 
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4.  Pressure  Relationships  for  Packed  Particle  Bed  (also  Cold  Frit): 


(Ergun  Relationship) 

f  150pvo(l  -  £)2 )  f  1 .75pve  |  v  |  (1  -  £) 

\  dp^  J  n  Dpe 


or: 


150jl(l  -  £)2 
D2pt?pA0W0 


1.75(1 -£) 
\  Dp?pA20 


1  f  A.+iO _ if  2 

.  p2^2  v  2AtA2  J  p,A,  ^  2AjA2 
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APPENDIX  C: 


MODEL  VALIDATION 

C.l  STEADY  STATE 

C.1.1  General 

In  order  to  provide  one  set  of  validations  for  the  simple  model,  it  may  be  compared 
to  the  reference  standard.  Steady  state  and  transient  results  for  both  models  are  compared 
graphically  in  the  following  subsections.  Values  for  reference  standard  parameters  were 
drawn  directly  from  T-l  and  values  for  the  simple  model  were  obtained  from  the  steady 
state  and  transient  source  codes  attached  as  Appendices  C  and  D.  Also,  each  model  used 
the  Sandia  National  Laboratory  PIPE  experiment  fuel  element  as  the  basis  for  its  overall 
configuration. 

C.1.2  Mass  Flowrate 

At  a  steady  state  condition,  coolant  mass  flowrate  will  vary  as  a  function  of  the  par¬ 
ticle  bed  power  density.  That  is,  the  mass  flowrate  will  adjust  itself  according  to  the 
resistance  to  flow  produced  by  the  element  as  a  whole  which  will  be  a  function  of 
temperature  and  pressure.  Figure  C.l  shows  how  mass  flowrate  varies  as  a  function  of 
power  density.  Keep  in  mind  that  the  average  temperature  of  the  fuel  and  the  exit  tem¬ 
perature  increase  with  an  increase  in  power  density.  Also  shown  in  Figure  C.l  are  the 
mass  flowrates  calculated  by  the  reference  standard  at  0,  1 ,  and  2  GW/m3.  There  is 
surprisingly  little  difference  between  the  simple  model  and  the  reference  standard  at  these 
power  densities. 
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Also  note  that  the  exit  plenum  0  is  adjusted  from  2.66  to  1.10  as  discussed  in  sec¬ 
tion  3.62.  This  gives  a  change  in  calculated  flow  (at  2  MW/m3)  from  3 1  g/s  (without  the 
adjustment)  to  37  g/s  (=15%  difference).  The  adjustment,  however,  did  enable  the 
behavior  of  the  whole  curve  to  be  well  represented. 

At  this  point,  it  must  be  noted  that  the  savings  in  computer  time  is  substantial.  The 
time  required  for  the  reference  standard  to  obtain  the  three  data  points  shown  on  Figure 
C.l  is  approximately  24  hours  (8  hours  per  point)  while  the  time  required  to  generate  the 
same  data  using  the  simple  model  is  approximately  12  seconds.  Further,  the  simple 
model  is  able  to  provide  many  more  data  points  in  the  region  of  interest. 
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C.13  Power 


The  rated  output  of  the  fuel  element  may  be  measured  as  the  thermal  power  of  the 
fuel  element.  This  is  calculated  by  multiplying  the  mass  flowrate  through  the  element  by 
the  change  in  enthalpy  across  the  fuel  element: 

Power  =  W  AH  C.l 

Figure  C.2  shows  the  relationship  between  fuel  element  thermal  power  to  power 
density  and,  as  expected,  the  relationship  is  linear  in  nature  and  expresses  a  heat  balance. 
The  calculated  value  of  thermal  power  at  2  MW/m3  is  plotted  for  the  reference  model  and 
corresponds  to  the  value  calculated  by  the  simple  model. 
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C.1.4  Exit  Temperature 

Figure  C.3  shows  how  the  temperature  of  the  coolant  exiting  the  fuel  element  varies 
with  power  density.  Because  the  exit  plenum  is  modeled  as  a  single  control  volume,  the 
exit  temperature  represents  an  average  temperature  of  the  coolant  exiting  the  element  and 
does  not  take  into  account  any  temperature  differences  which  may  occur  axially.  The  ref¬ 
erence  standard  has  five  axial  subdivisions  in  the  outlet  plenum,  each  of  which  may  have 
a  different  coolant  temperature.  Because  of  this  axial  separation  and  the  ability  of  the 
reference  standard  to  designate  different  axial  fuel  distributions,  the  temperatures  calcu¬ 
lated  in  the  exit  plenum  for  the  reference  standard  are  higher  at  the  lower  axial  position 
and  decreases  as  the  coolant  exits  the  plenum.  The  simple  model  can  not  make  this  same 
representation  and,  as  such,  exit  temperatures  may  be  slightly  different  between  the  two 
models  even  though  mass  flowrates,  power  output  and  enthalpy  changes  may  be  the 
same.  In  this  case,  values  for  mass  flowrate  at  various  power  density  levels  are  near  ref¬ 
erence  standard  values  but  are  not  exact. 

The  reference  temperature  plotted  in  Figure  C.3  is  the  average  outlet  plenum  tem¬ 
perature  at  2  and  2.1  MW/m3.  This  represents  a  difference  of  approximately  2.5%  and 
should  be  considered  acceptable  for  initial  fuel  element  design  but  final  analysis  should 
be  made  by  the  more  detailed  model. 

Velocity  is  measured  at  the  exit  of  the  plenum.  Because  exit  velocity  of  the  coolant 
is  inversely  proportional  to  the  density  of  the  coolant,  and  hence  direcdy  proportional  to 
the  temperature  of  the  coolant,  it  is  not  surprising  to  see  the  reference  standard  exit  veloc¬ 
ity  fall  below  the  curve  generated  by  the  simple  model.  Exit  velocity  is  shown  as  a  func¬ 
tion  of  power  density  in  Figure  C.4. 
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Steady  State  Exit  Temperature 
vs  Power  Density  for  Baseline 
PIPE  Fuel  Element 


Steady  State  Exit  Velocity 
vs  Power  Density  for  Baseline 
PIPE  Fuel  Element 


Velocity  (m/s) 


Simple  Model  Ref  Standard 


Figure  C.4:  Steady  State  Exit  Velocity  vs  Power  Density 
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C.1.5  Pressure 

The  calculations  of  sections  C.l  are  based  on  constant  inlet  pressure  and  constant 
outlet  pressure.  As  power  increases,  the  total  flow  through  the  fuel  element  drops  and  the 
fraction  of  the  total  pressure  drop  increases  in  the  control  volumes  at  the  higher  tempera¬ 
tures  Gower  densities).  C.5  shows  these  fraction  by  giving  the  calculated  pressures  at 
points  located  at  the  inlet  and  outlet  of  the  element  and  at  the  inlet  and  outlet  of  the  fuel. 
This  also  represents  the  boundaries  of  the  three  control  volumes. 

As  mentioned  in  T-l,  the  inlet  plenum  and  the  cold  frit  dominate  the  resistance  at  0 
GW/m3  but  as  power  density  (and  exit  temperature)  is  increased,  the  outlet  plenum 
assumes  a  greater  share  of  the  resistance.  The  resulting  pressure  drops  across  each  con¬ 
trol  volume  are  shown  in  Figure  C.6.  When  totaled,  the  sum  of  the  pressure  drops  will 
remain  a  constant  85  MPa,  which  is  established  as  an  initial  condition.  Figure  C.6  also 
shows  values  for  pressure  drop  across  the  outlet  plenum  in  the  reference  standard. 
Although  values  at  0  and  1  GW/m3  are  near  to  the  values  produced  by  the  simple  model, 
the  value  at  2  MW/m3  seems  to  be  inconsistent.  This  inconsistency  is  unexplained  and  is 
judged  to  be  unimportant  when  the  consistent  behavior  of  other  variables  is  noted. 
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Steady  State  Pressure  Drops 
vs  Power  Density  for  Baseline 
PIPE  Fuel  Element 


Pressure  (kPa) 


00  0.50  1.00  1.50  2.00 


Power  Density  (GW/m3) 


In  Plenum  &  CF 


Particle  Bed 


[ 


Out  Plenum  &  HF  fSM) 


Outlet  Plenum  (RF) 


Figure  C.6:  Steady  State  Pressure  Drops  vs  Power  Density 


C.2  NULL  TRANSIENTS 


As  a  test  for  validating  the  simple  model,  null  transients  were  analyzed  at  0, 1 ,  and 
2  MW/m3  and  subsequent  mass  flowrate  responses  were  examined.  Figures  C.7,  C.8  and 
C.9  show  that  there  is  a  smooth  behavior  throughout  the  transient  and  that  there  is  negli¬ 
gible  difference  between  the  mass  flowrate  before  the  transient  and  the  mass  flowrate 
after  the  transient.  Each  of  the  transients  start  at  .01  seconds  and  the  scale  for  flowrate  on 
the  y-axis  is  expanded  for  better  viewing. 
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Mass  Flowrate  Response  to  Null 
Transient  at  0  GW/m3  for 
Baseline  PIPE  Fuel  Element 

Flowrate  (g/s) 


0.01  0.012  0.014  0.016  0.018  0.02 


Time  (s) 

j  *  Transient  *  Steady  State 

i _ ! 

Figure  C.7:  Mass  Flowrate  Response  to  Null  Transient  at  0  GW/m3 
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C.3  TRANSIENTS 


In  Figures  C.10,  C.ll,  and  C.12,  calculations  for  a  1  second  transient  from  0 
GW/m3  to  2  GW/m3.  Calculations  for  the  simple  model  for  mass  flowrate,  thermal 
power,  and  exit  temperature  are  compared  to  the  reference  standard.  As  can  be  seen,  the 
reference  standard  reaches  the  final  set  of  values  quicker  than  the  simple  model.  Adjust¬ 
ment  could  be  made  to  correct  the  sluggishness  of  the  simple  model  (e.g.  by  increasing 
the  fuel-to-coolant  heat  transfer  parameter  and/or  by  artificially  decreasing  the  mass  of 
coolant  in  the  fuel  element).  These  adjustments  have  not  been  attempted  in  this  thesis 
study.  Both  the  simple  model  and  the  reference  standard  respond  quickly;  both  end  at 
essentially  the  same  condition. 
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Mass  Flowrate  Response  for  Is 
Transient  from  0  to  GW/m3  for 
Baseline  PIPE  Fuel  Element 


Flowrate  (g/s) 


Time  (s) 


Ref  Standard 


Simple  Model 


Transient  Starts  at  t=ls 


Figure  C.10:  Mass  Flowrate  Response  to  Baseline  1  s  Transient 


Thermal  Power  Ouput  Power  for  Is 
Transient  from  0  to  2  GW/m3  for 
Baseline  PIPE  Fuel  Element 


Exit  Temperature  Response  for  Is 
Transient  from  0  to  2  GW/m3  for 
Baseline  PIPE  Fuel  Element 


Temperature  (K) 
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Ref  Standard 


Transient  Starts  at  t=ls 


Figure  C.12:  Exit  Temperature  Response  to  Baseline  1  s  Transient 


C.4  NUMERICAL  STABILITY 

Accuracy  of  the  solution  of  the  control  volume  balance  equations  is  illustrated  in 
Figure  C.13.  That  is,  consider  that  the  balance  equations  are  written  in  terms  of  ordinary 
differential  equations  in  time.  Then  consider  solving  these  equations  (as  in  the  simple 
model)  by  discretizing  the  equations  in  the  time  variable.  Figure  C.13  gives  the  results 
related  to  the  errors  generated  in  this  discretization.  It  does  not  address  the  question  of 
discretization  in  space,  however. 

The  balance  equations  were  discretized  using  timesteps  of  100  ms,  10  ms,  and  1  ms. 
All  of  the  curves  are  very  close  and  do  not  show  much  variation.  The  plot  of  mass  flow- 
rate  corresponding  to  timesteps  of  1  ms  and  10  ms  seem  to  represent  the  solution  of  the 
differential  equations  better  than  the  curve  corresponding  to  100  ms.  The  results  for  100 
ms  indicate  an  acceptable  error  (less  than  5%)  and  could  be  used  for  faster  calculations. 
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APPENDIX  D: 


SOURCE  CODE  AND  DATA  ENTRY 

D.l  SOURCE  CODE  DESIGN 

The  source  code  for  the  analysis  of  the  response  of  a  particle  bed  reactor  fuel  ele¬ 
ment  is  provided  in  two  parts:  STEADY  and  UPPOWER.  STEADY.C  provides  a  means 
of  calculating  values  at  steady  state  conditions  for  a  variety  of  fuel  element 
configurations  and  operating  parameters.  This  enables  the  user  to  analyze  steady  state 
conditions  without  having  to  run  the  transient  program.  UPPOWER.C  is  similar  to 
STEADY.C  but  calculates  transient  values  over  a  specified  simulation  period. 

Both  source  codes  are  written  in  the  C  programming  language  and  compiled  using 
Microsoft™  Quick-C.  The  programs  are  intended  for  use  on  a  desktop  personal  computer 
equipped  with  a  hard  drive  but  may  run  on  a  computer  equipped  with  floppy  disk  drives 
if  one  drive  is  designated  as  the  "c  drive".  STEADY  produces  the  output  data  file 
STEADYST.DAT  and  UPPOWER  produces  the  data  file  POWER.DAT.  The  data  files 
generated  by  STEADY  and  UPPOWER  are  unformatted  files  which  are  easily  imported 
into  a  spreadsheet  (ie,  Lotus  123™)  which  then  allows  the  user  to  review  the  output 
graphically. 

D.2  DATA  ENTRY 

STEADY  and  UPPOWER  provide  fuel  element  geometry  and  initial  conditions 
which  are  similar  to  the  PIPE  experiments  as  initial  values  of  all  program  variables. 

These  values  may  be  changed  via  a  series  of  input  screens  presented  to  the  user  prior  to 


105 


beginning  an  analysis.  The  first  screen  (Figure  E.l)  initializes  all  important  geometry 
variables  which  may  be  changed  by  responding  to  the  prompts  under  the  list  of  parame¬ 
ters. 
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AFE  i  AR  AMETEFS  FOR  FPBR  FUEL  ELEMENT: 


THE  FOIL  -WING 


1 . 

Outer  Radius 

= 

0. 5445 SC 

VTi 

2  . 

Cold  Frit  Radius 

* 

0 .  C2?3*>0 

in 

3. 

Cold  Frit  Thickness 

* 

0.001*70 

TO 

4  . 

Dia  of  Cold  Frit  Particle 

- 

w.  00000270  in 

5  . 

Cold  Frit  Porosity 

= 

0 . 685300 

t . 

Particle  Bed  Thickness 

- 

0.012432 

m 

. 

Fuel  Particle  Diameter 

0 . 00C500 

m 

e. 

Particle  Bed  Void  Fraction 

- 

0.40000C 

9. 

Hot  Frit  Thickness 

* 

0 . 000760 

m 

10. 

Hot  Frit  Porosity 

- 

C. 230000 

n. 

Inlet  Manifold  Factor 

* 

0.950000 

12. 

Exit  Manifold  Factor 

* 

C. Ra5300 

13. 

Fuel  Element  Length 

- 

0.265000 

“ 

CHANGE  A  PARAMETER?  <1«Y  or  0-N/  2 

PLEASE  ENTER  THE  FAPAMETEP  NUMBER  AND  VALUE  (eg  12,  1.705): 

12,  1.0 

CHANGE  ANOTHER  PARAMETER  FOP  ELEMENT  A?  (2-Y  or  0-N)  0 

V _ 


Figure  E.l:  Data  Input  Screen  #1  for  STEADY 


Once  the  user  is  satisfied  with  the  fuel  element  geometry,  a  listing  of  these  values  as 
they  compare  to  the  baseline  PIPE  fuel  element  is  then  shown  on  the  screen  (Figure  E.2). 
This  is  followed  by  a  screen  which  allows  changes  to  specific  fuel  particle  properties 
(Figure  E.3)  and  a  final  input  display  which  allows  changes  to  general  operating  parame¬ 
ters  (Figure  E.4). 
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PARAMETER 

ELEMENT  A 

BASELINE 

1. 

Outer  Radius 

0.044580 

0 . 044580 

m 

2. 

Cold  Frit  Radius 

0.029360 

0.029360 

in 

3. 

Cold  Frit  Thickness 

0 . 001870 

0.001870 

m 

4  . 

Dia  of  Cold  Frit  Particle 

0 . 00000270 

0.00000270  m 

5  . 

Cold  Frit  Porosity 

0 . 685000 

0 . 685000 

6. 

Particle  Bed  thickness 

0.012432 

0.012432 

m 

7  . 

Fuel  Particle  Diameter 

0.000500 

0.000500 

m 

8  . 

Particle  Bed  Void  Fraction 

0.400000 

0 .400000 

9. 

Hot  Frit  Tv  .ckness 

0.000760 

0 . 000760 

m 

10 

.Hot  Frit  Porosity 

0.230000 

0.230000 

11 

.Inlet  Manifold  Factor 

0 . 950000 

0.950000 

12 

.Exit  Manifold  Factor 

0 . 888800 

0.888800 

13 

.Length  of  Element 

1.000000  m 

0.265000 

HIT  ANY  KEY  TO  CONTINUE 


Figure  E.2:  Data  Input  Screen  #2  for  STEADY 


r 

l . 

Fuel  Material 

OC2 

2 . 

Layer  1  Material 

= 

Low  D  Car 

3. 

Layer  2  Material 

High  D  CarZrC 

4  . 

Layer  3  Material 

ZrC 

5  . 

Radius  of  Fuel 

0.000117  m 

6  . 

Radius  of  Layer  1 

0.000150  m 

1 . 

Radius  of  Layer  2 

0.000200  in 

8  . 

Radius  of  Layer  3 

0.000250  m 

9  . 

Density  of  Fuel 

10500.000000  kg/m3 

10 

Density  of  Layer  1 

1000.000000  kg/m3 

11 

Density  of  Layer  2 

1900.000000  kg/m3 

12 

Density  of  Layer  3 

6300.000000  kg/m3 

13 

Cp  of  Fuel 

200.000000  J/kgK 

14 

Cp  of  Layer  1 

3000.000000  J/kgK 

15 

Cp  of  Layer  2 

3000.000000  J/kgK 

16 

Cp  of  Layer  3 

00.000000  J/kgK 

17 

k  for  Fuel 

** 

30.000000  W/m2K 

18 

k  for  Layer  1 

- 

1.500000  W/m2K 

19 

k  for  Layer  2 

- 

3.000000  W/m2K 

20 

k  for  Layer  3 

- 

40.000000  W/m2K 

CHANGE  A  PARAMETER?  ll-Y  or  0*N)0 


.'KAN  3£  AIM 7HEP  F Af  AMETEP 2  (i  =  Y  or  C'-li  i  v 


Figure  E.3:  Data  Input  Screen  #3  for  STEADY 
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THE  FOLLOWING  INITIAL  AND  FINAL  CONDITIONS  ARE  SET 


1. 

2. 

3. 

4. 

5. 

6. 
7. 


Initial  Inlet  Temperature  *  300.000000  K 
Initial  Inlet  Pressure  —  2000.000000  KPa 
Initial  Outlet  Pressure  =  1915.000000  KPa 
Initial  Power  Density  =  0.000000  GW/m3 
Final  Power  Density  -  2.500000  GW/m3 
Power  Density  Increment  =  0.100000  sec 
Initial  Guess  at  Mass  Flowrate  =  0.100000  kg/sec 


CHANGE  A  PARAMETER?  (1=Y  or  0-N)  1 


ENTER  PARAMETER  NUMBER,  VALUE  (eg  20,  40) 
3, 1900 


CHANGE  ANOTHER  PARAMETER?  (1-Y  or  0=N) 0 


Figure  E.4:  Data  Input  Screen  #4  for  STEADY 


UPPOWER  uses  the  same  data  input  screens  with  the  exception  of  input  screen  #4. 
UPPOWER  allows  the  user  to  select  a  value  for  the  initial  and  final  value  of  inlet  temper¬ 
ature,  inlet  pressure,  oudet  pressure  and  power  density  (Figure  E.5).  Further,  the  user 
may  choose  the  duration  of  each  transient.  All  transients  start  at  the  same  time  except  for 
inlet  temperature  which  has  a  1  s  delay  to  artificially  simulate  a  preheating  effect  from  a 
nuclear  thermal  rocket  nozzle  cooling  jacket. 


109 


THE  FOLLOWING  INITIAL  AND  FINAL  CONDITIONS  ARE  SET 


1 .  Initial  Inlet  Temperature 

2.  Final  Inlet  Temperature 

3.  Duration  of  Temperature  Change 

4 .  Initial  Inlet  Pressure 

5.  Final  Inlet  Pressure 

6.  Duration  of  Pressure  Change 

7.  Initial  Outlet  Pressure 

8.  Final  Outlet  Pressure 

9.  Duration  of  Pressure  Change 

10.  Initial  Power  Density 

11.  Final  Power  Density 

12.  Duration  of  Power  Density  Change 

13.  Time  Delay  to  Transient 

14.  Time  After  Transient 

15.  Time  Step 

16.  Initial  Guess  at  Mass  Flowrate 


300.000000  K 
300.000000  K 
1.000000  sec 
2000.000000  KPa 
2000.000000  KPa 
1.000000  sec 
1915.000000  KPa 
1915.000000  KPa 
1 . 000000  sec 
0.000000  GW/ m3 
2.000000  GW/ m3 
1 . 000000  sec 
1 . 000000  sec 
8.000000  sec 
0.010000  sec 
0.100000  kg/ sec 


CHANGE  A  PARAMETER?  (1-V  or  0-N) 0 

CHANGE  ANOTHER  PARAMETER?  (1-Y  or  0-N)0 
Data  Sample  Frequency? (print  every  xt h  data 
lO 


point . . ) 


Figure  E.5:  Data  Input  Screen  #4  for  UPPOWER 


It  should  be  noted  that,  as  a  final  input,  UPPOWER  asks  for  a  sampling  frequency. 
This  determines  the  amount  of  data  written  to  the  data  file.  Calculations  are  still  per¬ 
formed  at  the  required  timestep,  however.  For  example,  if  every  data  point  were  to  be 
written  to  the  data  file  using  a  timestep  of  1  ms  for  a  simulated  transient  lasting  10  s, 
there  would  be  10,000  data  points  written  to  the  output  file  (one  for  each  millisecond).  If 
a  sample  frequency  of  100  were  to  be  used  for  the  same  transient,  only  100  data  points 
would  be  written  to  the  output  file  (one  every  1 00  ms).  This  serves  only  to  reduce  the 
data  file  to  a  manageable  size  and  will  not  eliminate  any  of  the  calculations. 

D.3  Hydrogen  Equations  of  State 

To  analyze  fuel  element  transients,  properties  of  hydrogen  coolant  (such  as  density, 
viscosity  and  enthalpy)  must  be  determined  for  various  temperatures  and  pressures. 
Relationships  to  do  this  are  placed  at  the  end  of  the  source  code  as  subroutines  and  are 
accessed  during  the  main  part  of  the  program  as  needed. 
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The  density  of  hydrogen  is  determined  by  using  a  known  value  of  hydrogen  at  300 
K  and  101.325  kPa  (1  atmosphere)  and  using  the  ideal  gas  relationship  to  determine  other 
values.  The  subroutine  in  the  source  codes  determines  the  value  of  coolant  density  as  a 


function  of  temperature  and  pressure. 

Viscosity,  enthalpy,  coolant  heat  conductivity  and  Prantdl  number,  on  the  other 
hand  are  determined  as  a  function  of  temperature  only  since  the  contribution  due  to  a 
change  in  pressure  is  negligible  (R-l).  To  illustrate  this,  the  enthalpy  change  from  200  K 
to  1200  K  at  a  pressure  of  101.325  kPa  is  16406.2  kJ/kg  and  at  a  pressure  of  2  MPa  is 
16478.7  kJ/kg.  The  difference  between  the  two  values  is  approximately  0.4%.  This  is 
considered  negligible  for  the  simple  model.  The  other  properties  change  in  a  similar 
fashion. 
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D.4  SOURCE  CODE  FOR  STEADY 


/*  STEADY.C  *1 

/*  Copywrite  WILLIAM  E.  CASEY,  MAY  1990.  */ 

/*  The  author  hereby  grants  to  MIT  and  to  the  US  Government  */ 
/*  permission  to  reproduce  and  distribute  copies  of  this  code.  */ 

/*  This  source  code  was  written  using  the  Microsoft  QUICK-C  */ 
/*  programming  environment.  */ 

/*  Set  up  include  files:  */ 

#include  <c:\msc\include\stdio.h> 

#include  <c:\msc\include\conio.h> 

#include  <c:\msc\include\graph.h> 

#include  <c:Vnsc\include\stdUb.h> 

#include  <c:\msc\include\math.h> 

#define  PI  3.14159 

/*  Define  a  structure  for  geometry  related  variables:  */ 

struct  fuelelem 

{ 

float  OUTRAD;  /*  Outer  radius  of  Inlet  Plenum  */ 
float  CFRAD;  /*  Outer  radius  of  Cold  Frit  */ 
float  CFTHK;  /*  Cold  Frit  Thickness  */ 

float  CFDp;  /*  Diameter  of  Particles  in  Cold  Frit  */ 
float  CFPOR;  /*  Cold  Frit  Porosity  */ 

float  PBEDTHK;  /*  Particle  Bed  Thickness  */ 

float  PBEDDp;  f*  Diameter  of  Fuel  Particles  */ 
float  PBEDVOID;  /*  Particle  Bed  Void  Fraction  */ 
float  HFTHK;  /*  Hot  Frit  Thickness  */ 

float  HFPOR;  /*  Hot  Frit  Porosity  */ 

float  MFIN;  /*  Manifold  Factor  for  Inlet  Plenum  */ 
float  MFOUT;  f*  Manifold  Factor  for  Outlet  Plenum  */ 
float  LENGTH;  I*  Overall  Length  of  Element  */ 

}; 

/*  Define  a  structure  for  fuel  particle  composition:  */ 

struct  fuelpart 

{ 

charFUELMAT[10];  /*  Material  used  for  Fuel  */ 

char  L1MAT[10];  /*  Material  for  Layer  1  */ 

charL2MAT[10];  /*  Material  for  Layer  2  */ 

charL3MAT[10];  /*  Material  for  Layer  3  */ 

float  FUELRAD;  /*  Radius  of  Fuel  Core  */ 

float  LI  RAD;  /*  Radius  of  Layer  1  *1 

float  L2RAD;  /*  Radius  of  Layer  2  */ 

float  L3R  AD;  /*  Radius  of  Layer  3  */ 

float  FUELDEN;  /*  Density  of  Fuel  */ 

float  L1DEN;  f*  Density  of  Layer  1  *1 

float  L2DEN;  /*  Density  of  Layer  2  *1 

float  L3DEN;  /*  Density  of  Layer  3  *1 

float  FUELCp;  /*  Cp  for  Fuel  */ 

float  LICp;  /*  Cp  for  Layer  1  */ 

float  L2Cp;  /*  Cp  for  Layer  2  */ 

float  L3Cp;  /*  Cp  for  Layer  3  */ 

float  FUELK;  /*  Heat  Transfer  Coef  for  Fuel  */ 
float  L1K;  /*  Heat  Transfer  Coef  for  Layer  1  */ 

float  L2K;  /*  Heat  Transfer  Coef  for  Layer  2  */ 

float  L3K;  /*  Heat  Transfer  Coef  for  Layer  3  *1 

}; 
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/*  Define  a  structure  for  initial  and  final  conditions:  */ 

struct  conditions 

{ 

float  INITTIN;  /*  Initial  Value  of  T  in  */ 

float  INIT_PIN;  f*  Initial  Value  of  P  in  */ 

float  INIT_POUT;  /*  Initial  Value  of  P  out  */ 

float  INIT_PRDEN;  /*  Initial  Value  of  Power  Density  */ 
float  FIN_PRDEN;  t*  Final  Value  of  Power  Density  */ 
float  delta_t;  /*  Increment  */ 

float  INIT_W;  /*  Initial  Guess  at  Mass  Flowrate  */ 

); 

void  showdefaults(  struct  fuelelem  *e_ptr ); 

void  fuelelemsumf  struct  fuelelem  *e_ptr,  struct  fuelelem  *f_ptr ); 

void  showfuelpart  (  struct  fuelpart  *g_ptr ); 

void  showconditions  (  struct  conditions  *h_ptr ); 

/*  Define  H2  state  relationships  found  at  end  of  source  code:  */ 

float  h2density  (float  xx,  float  yy); 

float  h2viscosity  (float  xx); 

float  h2heatxfer  (float  xx); 

float  h2cp  (float  xx); 

float  h2prandle  (float  xx); 

float  h2H  (float  xx); 

/*  **  BEGINNING  OF  COMPUTATIONAL  CODE  ************************  */ 
main() 

{ 

/*  Define  variables  for  use  in  source  code:  */ 

double  WO,  Wl; 
double  pinl,  pinll,  pinlll,  pout; 
double  pbarl,  pbarll,  pbarlll; 
double  Tin,  Tout,  Tbarll; 

float  Pinl,  Pinll,  Pinlll,  Pout,  Pbarl,  Pbarll,  Pbarlll; 

float  uin,  uout,  ubarll; 

float  RI,  Rfl,  Rm,  RTOTAL; 

float  PDROPII,  PDROPH1,  PDROPHU; 

float  PDROPI2,  PDROPII2,  PDROPD32; 

float  Hinll,  Hinlll,  HinflU,  Houtl; 

float  HinI2,  HinII2,  HinIII2,  Hout2; 

float  Hbarll,  HbarI2,  Hbarlll,  Hbarffi,  HbarlUl,  HbarIII2; 

float  Qbed; 

float  AinI,  AoutI,  ACF,  Ainll,  Aoutll,  AHF,  AinlH,  Aoutlll; 
float  Abarll; 

float  CPin,  CPout,  CPtemp; 
float  Werror,  CPerr, 
float  Deq_in,  Deq_out; 
float  Rif,  RIa,RIk,  Rim; 
float  RCF,  Rel,  Rem,  Rlla; 
float  Rfflf,  RHIa,  Rlllk,  RDIm; 
float  mfout; 

float  powerdensity,  vout; 
int  ch,  chr, 
float  aa,  alpha,  cc; 
float  Htemp,  Herror, 
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static  struct  fuelelem  A  = 

{  .04458,  .02936,  .00187,  .0000027,  .685,  .012432,  .0005, 

.4,  .00076,  .23,  .95, 1.1,  .265 

); 

static  struct  fuelelem  B  = 

{  .04458,  .02936,  .00187,  .0000027,  .685,  .012432,  .0005, 

.4,  .00076,  .23,  .95, 1.1,  .265 

}; 

static  struct  fuelpart  C  = 

{  "UC2",  "Low  D  Car”,  "High  D  Car",  "ZrC",  .000117,  .00015, \ 

.0002,  .00025,  10500, 1000,  1900,  6300, 200,  3000,  3000,  \ 

200,  30, 1.5,3,40 

); 

static  struct  conditions  D  = 

{  300,2000, 1915,  0,2.5,  .1,  .1 

); 

FILE  *out; 

out  =  fopen  (  ”c:steadyst.dat","w+"  ); 
ch  =  1; 

while  (ch  =  1  ) 

{ 

sfaowde faults  ( <fcA  )* 

printfOViCHANGE  ANOTHER  PARAMETER  FOR  ELEMENT  A?  (1=Y  or  0=N)"); 
scanf("%i",  &ch); 

} 

fuelelemsum  (  &A,  &B  ); 
ch=  1; 

while  (ch  =  1 ) 

{ 

showfuelpart  (  &C  ); 

printf(’NtiCHANGE  ANOTHER  PARAMETER?  (1=Y  or  0=N)"); 
scanf("%i",  &ch); 

) 

ch  =  1; 

while  (ch  ==  1  ) 

{ 

showconditions  ( &D ); 

printf('ViCHANGE  ANOTHER  PARAMETER?  (1=Y  or  0=N)"); 
scanf(”%i",  &ch); 

1 

/*  Determination  of  Steady  State  Values-Initial  Conditions  */ 

fprintf(  out, "STEADY  STATE  VALUES  FOR  THE  FOLLOWING  INITIAL  CONDITlONSSn"); 
fprintf(  out,"\nInlet  Temp  =  %.2f  KNn",  D.INIT_TIN); 
fprintf(  out, "Inlet  Press  =  %.2f  kPaNn",  D.INIT_PIN); 
fprintf(  out,"Outlet  Press  =  %.2f  kPaNnNn”,  D.INIT.POUT); 

fprintf(out,”PWRDN\tQbe<NW\tPDI\tPDir^PDniMToutNtTbarn\tPmI\tPinIl\tPinIIINtPoutMvoutVi\n"); 
W1  =  D.INTT_W; 

W0  =  Wl; 

AinI  =  PI  *  (A.OUTRAD  *  A.OUTRAD  -  A.CFRAD  *  A.CFRAD); 

AoutI  =  2  *  PI  *  A.CFRAD  *  A.LENGTH; 

ACF  =  AoutI  *  A.CFPOR- 

AinE  =  2  *  PI  *  (A.CFRAD  -  A.CFTHK)  *  A.LENGTH  *  A.PBEDVOID; 

Aoufll  =  2  *  PI  *  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK)  *  A.LENG'fH  *  \ 

A.PBEDVOID; 

Abarll  =  .5  *  (Anil  +  AouUI); 
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Ainlll  =  2  *  PI  *  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK  -  A.HFTHK)  *  \ 

A.LENGTH; 

Aoutm  =  PI  *  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK  -  A.HFTHK)  *  \ 

(A.CFRAD  -  A.CFTHK  -  A.PBEDTHK  -  A.HFTHK); 

AHF  =  Ainlll  *  A.HFPOR; 

Deq_in  =  4  *  AinI  1(2  *  PI  *  (A.OUTRAD  +  A.CFRAD)); 

Deq_out  =  4  *  Aoutm  /  (2  *  PI  *  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK  -  \ 

A.HFTHK)); 

for  (powerdensity  =  D.INIT_PRDEN;  powerdensity  <  D.FIN_PRDEN+.1;\ 
powerdensity  +=  .1) 

{ 

Qbed  =  powerdensity  *  1000000000  *  PI  *  ((A.CFRAD  -  A.CFTHK)  *  \ 

(A.CFRAD  -  A.CFTHK)  -  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK)  *  \ 
(A.CFRAD  -  A.CFTHK  -  A.PBEDTHK))  *  A.LENGTH; 

Tin  =  D.INIT_TIN; 

Hinll  =  h2H(Tin); 

PinI  =  D.INIT_PIN; 

Pinll  =  D.INIT_PIN; 

Pinffl  =  D.INIT_POUT ; 

Pout  =  D.INIT_POUT; 

Tout  =  Tin; 
do 
{ 

W0  =  Wl; 

Houtl  =  Q  bed/WO  +  Hinll; 
do 
{ 

Htemp  =  h2H(Tout); 

Herror  =  Houtl  -  Htemp; 

Tout  =  Tout  +  Henor/100000; 

}  while  ( Herror/lOOOOOO  >  .000001  ); 

/*  Set  up  all  variables  */ 

Tbarll  =  ,5*(Tin  +  Tout); 

Pbarl  =  .5*(PinI  +  Pinll); 

Pbarll  =  ,5*(PinII  +  Pinm); 

Pbarin  =  ,5*(Pinm  +  Pout); 
pbarl  =  h2density(Tin,  Pbarl); 
pbarll  =  h2density  (Tbarll,  Pbarll); 
pbarIII=  h2density(Tout,  Pbarlll); 
uin  =  h2viscosity(Tin); 
ubarll  =  h2viscosity(TbaiH); 
uout  =  h2viscosity(Tout); 
pinll  =  h2density(Tin,  Pinll); 
pinlll  =  h2density(Tout,  PinHI); 
pout  =  h2density(Tout,Pout); 

/*  Inlet  Resistance  Calculations  */ 

Rel  =  W0*Deq_in  /  (uin*  AinI); 

Rif  =  ,138*pow(ReI,  -,151)*(A.LENGTH/Deq_in)/(2*pbarI*  AinI*  AinI); 

RIa  =  (l/(AinII*AinII)  -  l/(AinI*AinI))/(2*pbarI); 

RIk  =  (l-(AinI/AoutI))*(l-(AinI/AoutI))/(2*pbarI*AinI*AinI)  \ 

+  .5  *  (l-(ACF/AoutI))*(l-(ACF/AoutI))/(2*pinn*ACF*ACF)  \ 

+  .5  *  ( l-( AinII/ACF))*( l-(AinH/ACF))  /  (2*pinII*Ainn*Ainn); 

Rim  =  A.MFIN  /  (pbari*  AinI*  AinI); 

Rim  =  Rim  +  A.MFIN^pbar^Aoutl*  AoutI); 

RCF  =  150*uin*(l-A.CFPOR)*(l-A.CFPOR)/(A.CFDp*A.CFDp*pinn*Aoud*W0*  \ 
A.CFPOR*  A.CFPOR*  A.CFPOR)  \ 
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+  1.75*(l-A.CFPOR)/(A.CFDp*pinD*AoutI*AoutI*A.CFPOR*A.CFPOR*  \ 
A.CFPOR); 

RCF  =  RCF  *  A.CFTHK; 

RI  =  RIf/2  +  RIa  +  RDc  +  Rim  +  RCF; 

/*  Resistance  Calculations  for  the  Particle  Bed  */ 

RII  =  150*ubarII*(l-A.PBEDVOID)*(l-A.PBEDVOID)/(A.PBEDDp*A.PBEDDp*\ 
pbarn*Abarn*WO*A.PBEDVOID*A.PBEDVOID)  \ 

+  1.75*(l-A.PBEDVOID)/(A.PBEDDp*pbarn*Abarn*AbarII*A.PBEDVOID); 
RH  =  RII  *  A.PBEDTHK; 

Rlla  =  ( l/(pinm*AinIII)  -  l/(pinn*Ainn)  )*((Ainn+Ainm>/  \ 

(2*Ainn*Ainin)); 

RU  =  RII  +  RHa; 

/*  Resistance  for  Exit  Plenum  */ 

Rein  =  W0*Deq_out  /  (uout*AoutIII); 

RTHf  =  .138*pow(Rein,  -.151)*(A.LENGTH/Deq_out)/(2*pbarIII*AouaiI\ 

*Aoutffl); 

Rllla  =  (l/(Aoutm*Aoutm)-l/(Aoutn*Aoutn))/(2*pbarm); 

Rink  =  .5*(l-AHF/Aoutn)*(l-AHF/Aoutn)/(2*pinm*AHF*AHF)\ 

+  (l-AHF/Ainm)*(l-AHF/Ainm)/(2*pinm*AHF*AHF)  \ 

+  (l-Aoutm/Ainm)*(l-Aoutm/Ainm)/(2*pbarm*  \ 

AoutIU*Aoutni); 
mfout  =  A.MFOUT; 

Rlllm  =  mfout  /  (pbarin*Aoutm*Aoutin)  +  mfout  /  (pbarHI  \ 

*AinIH*Ainin); 

RHI  =  Rnif/2  +  RIHa  +  RMk  +  REIm; 

/*  Final  Assembly  */ 

RTOTAL  =  RI  +  Rn  +  RHI; 

W1  =  (PinI  -  Pout)*  lOOO/RTOTAL; 

W1  =  sqit(Wl); 

Werror  =  WO- Wl; 

Werror  =  fabs(Werror); 

PDROPI1  =  RI*W1*W1/1000; 

PDROPH1  =  Rn*Wl*Wl/1000; 

PDROPmi  =  Rffl*WI*Wl/1000; 

PinH  =  Pinl-PDROPIl; 

PinDI  =  Pout+PDROPm  1 ; 
vout  =  Wl/(pout*Aoutni); 

}  while  (Werror  >  .0000001); 

W1=W1*1000; 

Qbed  =  Qbed/1000; 

fprintf(out,"%.2N%.lfNt%.5Pt%.2N%.2Pt%.2fM%.2fNt%.2N%.2fNt%.2N%.2fM%.2fM%.2Po",\ 
powerdensity,  Qbed,  W 1  ,PDROPI  1  ^DROPH 1  .PDROPHI 1  ,Tout,TbarII,Pi- 
nIJPinn,PinIII,Pout,vout); 

printf  ("Pwrden  =  %fStW  =  %fSn",  powerdensity,  Wl); 

Qbed  =  Qbed*  1000; 

W1=W1/1000; 
if  (Tout  >2700) 

goto  meltdown; 

1 

meltdown: 
if  (Tout>2700) 

{  printfC’Tout  >  2700"); 

fprintf(out,"Temp  exceeds  2700"); 
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/****>,***********************************»***************************/ 
void  showdefaults(  struct  fuelelem  *e_ptr ) 

{ 

int  num; 
float  entry; 
int  c; 

clearscreen(O)' 

printf("THE  FOLLOWING  ARE  PARAMETERS  FOR  PPBR  FUEL  ELEMENT :\n\n”); 
printf("  1 .  Outer  Radius  =  %f  m\n",  e_ptr->OUTRAD); 

printf("2.  Cold  Frit  Radius  =  %f  m\ri',  e_ptr->CFRAD); 

printf("3.  Cold  Frit  Thickness  =  %f  m\n",  e_ptr->CFTHK); 
printf("4.  Dia  of  Cold  Frit  Particle  =  %.8f  mV',  e_ptr->CFDp); 
printf("5.  Cold  Frit  Porosity  =  %fc",  e_ptr->CFPOR); 
printf("6.  Particle  Bed  Thickness  =  %f  m\n",  e_ptr->PBEDTHK); 
printf(  "7.  Fuel  Particle  Diameter  =  %{  m\n”,  e_ptr->PBEDDp); 
printf(”8.  Particle  Bed  Void  Fraction  =  %fSri\  e_ptr->PBEDVOID); 
printf("9.  Hot  Frit  Thickness  =  %f  m\n",  e_ptr->HFTHK); 
printf("10.  Hot  Frit  Porosity  =  %fW,  e_ptr->HFPOR); 
printf("  1 1 .  Inlet  Manifold  Factor  =  %fSn”,  e_ptr->MFIN); 

printf("12.  Exit  Manifold  Factor  =  e_ptr->MFOUT); 
printf("13.  Fuel  Element  Length  =  %f  mVi",  e_ptr->LENGTH); 
printf('%nVi\nCH AN GE  A  PARAMETER?  (1=Y  or  0=N)"); 
scanf(  "%i",  &c ); 
if  (  c  —  0  ) 
return; 

printf("\nPLEASE  ENTER  THE  PARAMETER  NUMBER  AND  VALUE  (eg  12, 1.705):W); 
scanf(  "%i,  %f &num,  &entry  ); 
switch  (num) 

{ 

case  1: 

e_ptr->OUTRAD  =  entry; 
break; 

case  2: 

e_ptr->CFRAD  =  entry; 
break; 

case  3: 

e_ptr->CFTHK  =  entry; 
break; 
case  4: 

e_ptr->CFDp  =  entry; 
break; 

case  5: 

e_ptr->CFPOR  =  entry; 
break; 

case  6: 

e_ptr->PBEDTHK  =  entry; 
break; 

case  7: 

e_ptr->PBEDDp  =  entry; 
break; 

case  8: 

e_ptr->PBEDVOID  =  entry; 
break; 

case  9: 

e_ptr->HFTHK  =  entry; 
break; 
case  10: 
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e_ptr->HFPOR  =  entry; 
break; 
case  11: 

e_ptr->MFIN  =  entry; 
break; 
case  12: 

e_ptr->MFOUT  =  entry; 
break; 
case  13: 

e_ptr->LEN  GTH  =  entry; 
break; 
default: 

printf('\n\nTRY  ANOTHER  PARAMETER^"); 

_clearscreen(0); 

break; 

} 

return; 

} 

/*  ******************************************************************  * j 

void  fuelelemsum  ( struct  fueMem  *e_ptr,  struct  fuelelem  *f_ptr ) 

_clearscreen  (0); 

printf  ("PARAMETERNM'tELEMENT  AM\tBASELINE\ri'o"); 
printf("l.  Outer  Radius\N\t%fvt\t%f  m\n",  e_ptr->OUTRAD,  f_ptr->OUTRAD); 
printf("2.  Cold  Frit  RadiusNtNt%fNtNt%f  mVt",  e_ptr->CFRAD,  f_ptr->CFRAD); 
printf(”3.  Cold  Frit  Thickness\f\t%NNt%f  mNn",  e_ptr->CFrHK,  f_ptr->CFTHK); 
printf("4.  Dia  of  Cold  Frit  Particle\t%.8fNftt%.8f  m'n",  e_ptr->CFDp,  f_ptr->CFDp); 
printf(”5.  Cold  Frit  Porosity\lNt%fNtNt%f\n",  e_ptr->CFPOR,  f_ptr->CFPOR); 
printf("6.  Particle  Bed  thickness\t%fM\t%f  m\n",  e_ptr->PBEDTHK,  f_ptr->PBEDTHK); 
printf("7.  Fuel  Particle  Diameter\t%fNtNt%f  mNn”,  e_ptr->PBEDDp,  f_ptr->PBEDDp); 
printf("8.  Particle  Bed  Void  Fraction\t%NNt%f\n”,  e_ptr->PBEDVOID,  f_ptr->PBEDVOID); 
printf(”9.  Hot  Frit  ThicknessNtNt%NNt%f  mNn",  e_ptr->HFTHK,  f_ptr->HFTHK); 
printf("10.Hot  Frit  Porosity\fa%P*Nt%fNn",  e_ptr->HFPOR,  f_ptr->HFPOR); 
printf("  11. Inlet  Manifold  Factoi\t%f\lNt%f\n",  e_ptr->MFIN,  f_ptr->MFIN); 
printf('T2.Exit  Manifold  Factoi\t't%fNtNt%f\n",  e_ptr->MFOUT,  f_ptr->MFOUT); 
printf("13.Lengtb  of  Elemenl\t\t%f  mNtM%f*i",  e_ptr->LENGTH,  f_ptr->LEN  GTH); 
printf("NnNnHIT  ANY  KEY  TO  CONTINUE"); 
getchQ; 


} 

I*  ******************************************************************  */ 

void  showfuelpart  ( struct  fuelpart  *g  ptr ) 

{ 

int  num; 
float  entry; 
int  c; 

_clearscreen(0); 

printf('THE  FOLLOWING  ARE  PARAMETERS  FOR  THE  FUEL  PARTICLENn"); 

printf("l.  Fuel  Material  =  %sNn",  g_ptr->FUELMAT); 

printf("2.  Layer  1  Material  =  %sNn",  g_ptr->LlMAT); 

printf(”3.  Layer  2  Material  =  %sNn”,  g_ptr->L2MAT); 

printf("4.  Layer  3  Material  =  %sNn”,  g_ptr->L3MAT); 

printf("5.  Radius  of  Fuel  =  %f  mV,  g_ptr->FUELRAD); 

printf("6.  Radius  of  Layer  1  =  %f  mNn",  g_ptr->LlRAD); 

printf("7.  Radius  of  Layer  2  =  %f  mNn",  g_ptr->L2RAD); 

printff'8.  Radius  of  Layer  3  =  %f  mNn",  g_ptr->L3RAD); 

printf("9.  Density  of  Fuel  =  %f  kg/m3V,  g_ptr->FUELDEN); 

printf("10.Density  of  Layer  1  =  %f  kg/m3Nn",  g_ptr->LlDEN); 


118 


printf("  11. Density  of  Layer  2  =  %f  kg/m3\n”,  g_ptr->L2DEN); 
printf("12.Density  of  Layer  3  =  %f  kg/m3\n",  g_ptr->L3DEN); 
printf("13.  Cp  of  Fuel  =  %f  J/kgKNn",  g_ptr->FUELCp); 

printf("14.  Cp  of  Layer  1  =  %f  J/kglOn",  g_ptr->LlCp); 
printf("15.  Q>  of  Layer  2  =  %f  J/kgK\n",  g_ptr->L2Cp); 
printf("16.  Cp  of  Layer  3  =  %f  J/kgKNn",  g_ptr->L3Cp); 
nrintf(”17.  k  for  Fuel  =  %f  W/m2K\n",  g_ptr->FUELK); 
printf("18.  k  for  Layer  1  =  %f  W/m2K>o",  g_ptr->LlK); 
printf("19.  k  for  Layer  2  =  %f  W/m2K>o”,  g_ptr->L2K); 
printf(''20.  k  for  Layer  3  =  %i  W/m2KSn",  g_ptr->L3K); 
printf(’\nCHANGE  A  PARAMETER?  (1=Y  or  0=N)"); 
scanf(  "%i",  &c ); 
if  ( c  ==  0  ) 


return; 

printf("\nENTER  PARAMETER  NUMBER  (  NOT  1-4  ),  VALUE  (eg  20, 40)W); 
scanf(  ”%i,  %€',  &num,  &entry  ); 
switch  (num) 

{ 

case  5: 

g_ptr->FUELRAD  =  entry; 
break; 

case  6: 

g_ptr->LlRAD  =  entry; 
break; 
case  7: 

g_ptr->L2RAD  =  entry; 
break; 

case  8: 

g  ptr->L3RAD  =  entry; 
break; 

case  9: 

g_ptr->FUELDEN  =  entry; 
break; 
case  IP- 

g_ptr->LlDEN  =  entry; 
break; 
case  11: 

g  ptr->L2DEN  =  entry; 
break; 
case  12: 

g  ptr->L3DEN  =  entry; 
break; 
case  13: 

g_ptr->FUELCp  =  entry; 
break; 
case  14: 

g_ptr->LlCp  =  entry; 
break; 
case  15: 

g_ptr->L2Cp  =  entry; 
break; 
case  16: 

g  ptr->L3Cp  =  entry; 
break; 
case  17: 

g  ptr->FUELK  =  entry; 
break; 
case  18: 


119 


g_ptr->LlK  =  entry; 
break; 
case  19: 

g_ptr->L2K  =  entry; 
break; 
case  20: 

g_ptr->L3K  =  entry; 
break; 
default: 

printfC’CANNOT  CHANGE  THIS  PARA  METER... SORRY.Nn"); 

getch(); 

break; 

) 

return; 

} 

I*  it:*********************************:********************************  *  I 

void  showconditions  (  struct  conditions  *h_ptr ) 

{ 

int  num; 
float  entry; 
int  c; 

_clearecreen(0); 

printf("THE  FOLLOWING  INITIAL  AND  FINAL  CONDITIONS  ARE  SEToW); 
printf("l.  Initial  Inlet  Temperature  =  %f  K\n”,  h_ptr- >INIT_TIN); 
printf(”2.  Initial  Inlet  Pressure  =  %f  KPavo",  h_ptr->INIT_PIN); 
printf("3.  Initial  Outlet  Pressure  =  %f  KPaNn",  h_ptr->INn'_POUT); 
printf("4.  Initial  Power  Density  =  %f  GW/m3\n",  h_ptr->INIT_PRDEN); 

printf("5.  Final  Power  Density  =  %f  GW/m3\n",  h_ptr->FIN_PRDEN); 

printf("6.  Power  Density  Increment  =  %f  secW,  h_ptr->delta_t); 
printf("7.  Initial  Guess  at  Mass  Flowrate  =  %f  kg/secV,  h_ptr->INTT_W); 
printf("SnCHANGE  A  PARAMETER?  (1=Y  or  0=N)”); 
scanf(  "%i",  &c  ); 
if  ( c  =====  0  ) 
return; 

printf("\nENTER  PARAMETER  NUMBER,  VALUE  (eg  20, 40)W); 
scanf(  "%i,  %f ',  &num,  &entry  ); 
switch  (num) 

{ 

case  1: 

h_ptr->INIT_TIN  =  entry; 
break; 

case  2: 

h_ptr->INIT_PIN  =  entry; 
break; 

case  3: 

h_ptr->INIT_POUT  =  entry; 
break; 
case  4: 

h_ptr->INIT_PRDEN  -  entry; 
break; 

case  5: 

h_ptr->FTN_PRDEN  =  entry; 
break; 

case  6: 

h_ptr->delta_t  =  entry; 
break; 

case  7: 

h_ptr->INIT_W  =  entry; 
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break; 

default: 

printf("CANNOT  CHANGE  THIS  PARAMETER.. .SORRY.Nn "); 

getch(); 

break; 

) 


return; 


float  h2density(float  xx,  float  yy) 

{ 

float  dens; 

dens  =  .08185  *  300  /  xx  *  yy  /  101.325; 
return  dens; 

} 

I*  ******************************************************************  * j 

float  h2viscosity(float  xx) 
float  vise; 

vise  =  .000008963  *  pow(  xx/300,  .6733); 
return  vise; 

} 

j*  ******************************************************************  *  j 

float  h2heatxfer( float  xx) 

{ 

float  teml; 
int  i; 

float  tem2; 

float  heat[45]  =  {0,  0362,  .0665,  .0981,  .1282,  .1561,  .182,  .206,  ,228,\ 

.251,  .272,  .292,  .315,  .333,  .351,  .3665,  .384,  .398,\ 

.412,  .426,  .44,  .452,  .464,  .476,  .488,  .500.  .512,  \ 

.524,  .536,  .548,  .560,  .572,  .584,  .5%,  .608,  .62,  \ 

.632,  .644,  .656,  .668,  .680,  .692,  .704,  .716,  .728  }; 

tem2  =  xx  /  50.0; 
i  =  floor(tem2); 

teml  =  heat[i]  +  (heat[i+l]  -  heat[i])  *  ((xx  -  (i  *  50))  /  50 ); 
return  teml; 

I 

ft*****************************************************************  +1 

float  h2cp(float  xx) 

{ 

float  temcp; 
if  ( xx  <=  420  ) 

temcp  =  4.1868  *  1000  *  (1.5395  +  .0150825  *  xx  -  4.02449e-5  *  \ 
xx  *  xx  +  3.63544e-8  *  xx  *  xx  *  xx); 

else 

temcp  =  4.1868  *  1000  *  (3.58927  -  5.55096e-4  *  xx  +  6.94235e-7  *  xx\ 
*  xx  -  1 .45 155e- 10  *  xx  *  xx  *  xx); 

return  temcp; 

} 

I*  ****************************************************************** 

float  h2H(float  xx) 

{ 

float  temH; 
float  temHl; 
temH=0; 
temH  1=0; 
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if  (xx  <=  420) 

I 

temH=4.1868*1000*(1.5395*xx  +  ,0150825*xx*xx/2  -  4.02449e-5*  \ 

xx*xx*xx/3  +  3.63544e-8*xx*xx*xx*xx/4); 
temH=temH  -  4. 1868*  1000*(  1.5395*50  +  .0150825*50*50/2  -  4.02449e-5*  \ 
50*50*50/3  +  3.63544e-8*50*50*50*50/4); 

) 

else 

{ 

temHl=4.1868*1000*(3.58927*xx  -  5.55096e-4*xx*xx/2  +  6.94235e-7*  \ 

xx*xx*xx/3  -  1.45155e-10*xx*xx*xx*xx/4); 
temHl=temHl  -  4.1868*1000*(3.58927*420  -  5.55096e-4*420*420/2  +  \ 

6.94235e-7*420*420*420/3  -  1. 45 155e- 10*420*420*420*420/4); 
temH=temHl  +  4.1868*1000*(1.5395*420+.0150825*420*420/2-4.02449e-5*  \ 
420*420*420/3  +  3.63544e-8*420*420*420*420/4); 
temH=temH  -  4.1868*1000*(1.5395*50  +  .0150825*50*50/2  -  4.02449e-5*  \ 
50*50*50/3  +  3.63544e-8*50*50*50*50/4); 

} 


return  temH; 

} 

j *  ******************************************************************  * j 

float  h2prandle(float  xx) 

{ 

float  teml; 
int  i; 

float  tem2; 

float  pran[23]  =  (.000,  .712,  .719,  .706,  .690,  .675,  .664,  \ 

.659,  .664,  .676,  .686,  .703,  .715,  .733,  \ 

.748,  .763,  .778,  .793,  .808,  .823,  .838,  \ 

.853,  .868); 

tem2  =  xx  / 100.0; 
i  =  floor(tem2); 

teml  =  pranfi]  +  (pran[i+l]  -  pran[i])  *  ((xx  -  (i  *  100))  /  100 ); 
return  teml; 

1 

j*  ****************************************************************** 
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E.5  SOURCE  CODE  FOR  UPPOWER 


/*  UPPOWER.C  */ 

/*  Copy  write  WILLIAM  E.  CASEY,  MAY  1990.  */ 

/*  The  author  hereby  grants  to  MIT  and  to  the  US  Government  */ 
/*  permission  to  reproduce  and  distribute  copies  of  this  code.  */ 

/*  This  source  code  was  written  using  the  Microsoft  QUICK-C  */ 
/*  programming  environment.  */ 

/*  Set  up  include  files:  */ 

#include  <c:\msc\include\stdio.h> 

#include  <c:Nmsc\include\conio.h> 

#include  <c:\msc\include\graph.h> 

#include  <c:\msc\include\stdlib.h> 

#include  <c:\msc\include\math.h> 

#define  PI  3.14159265442 
struct  fuelelem 

I 

float  OUTRAD;  /*  Outer  radius  of  Inlet  Plenum  */ 
float  CFRAD;  f*  Outer  radius  of  Cold  Frit  */ 
float  CFTHK;  /*  Cold  Frit  Thickness  */ 

float  CFDp;  /*  Diameter  of  Particles  in  Cold  Frit  */ 
float  CFPOR;  f*  Cold  Frit  Porosity  */ 

float  PBEDTHK;  /*  Particle  Bed  Thickness  */ 

float  PBEDDp;  /*  Diameter  of  Fuel  Particles  */ 
float  PBEDVOID;  /*  Particle  Bed  Void  Fraction  */ 
float  HFTHK;  /*  Hot  Frit  Thickness  */ 

float  HFPOR;  /*  Hot  Frit  Porosity  */ 

float  MFIN;  /*  Manifold  Factor  for  Inlet  Plenum  */ 
float  MFOUT;  /*  Manifold  Factor  for  Outlet  Plenum  */ 
float  LENGTH;  f*  Overall  Length  of  Element  */ 

}; 

struct  fuelpart 

{ 

char  FUELMAT  [10];  /*  Material  used  for  Fuel  */ 

char  L1MAT[10];  /*  Material  for  Layer  1  */ 

char  L2MAT[10];  /*  Material  for  Layer  2  */ 

char  L3MAT[10];  /*  Material  for  Layer  3  */ 

float  FUELRAD;  /*  Radius  of  Fuel  Core  */ 

float  L1R  AD;  /*  Radius  of  Layer  1  */ 

float  L2RAD;  /*  Radius  of  Layer  2  */ 

float  L3RAD;  /*  Radius  of  Layer  3  */ 

float  FUELDEN;  /*  Density  of  Fuel  */ 

float  LI  DEN;  /*  Density  of  Layer  1  */ 

float  L2DEN;  /*  Density  of  Layer  2  */ 

float  L3DEN;  /*  Density  of  Layer  3  */ 

float  FUELCp;  /*  Cp  for  Fuel  */ 

float  LI  Cp;  /*  Cp  for  Layer  1  */ 

float  L2Cp;  /*  Cp  for  Layer  2  */ 

float  L3Cp;  /*  Cp  for  Layer  3  */ 

float  FUELK;  /*  Heat  Transfer  Coef  for  Fuel  */ 
float  LI  K;  /*  Heat  Transfer  Coef  for  Layer  1  */ 

float  L2K;  /*  Heat  Transfer  Coef  for  Layer  2  */ 

float  L3K;  /*  Heat  Transfer  Coef  for  Layer  3  */ 

}; 

struct  conditions 

( 
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*/ 


float  INIT_TIN;  /*  Initial  Value  of  T  in  */ 

float  HN_TIN;  /*  Final  Value  of  T  in  */ 

float  DUR_TIN;  /*  Duration  of  T  in  Transient 
float  INIT_PIN;  /*  Initial  Value  of  P  in  */ 

float  FIN.PIN;  /*  Final  Value  of  P  in  */ 

float  DUR_PIN;  /*  Duration  of  P  in  Transient  */ 
float  INITJPOUT;  /*  Initial  Value  of  P  out  */ 

float  FIN_POUT;  /*  Final  Value  of  P  out  */ 

float  DUR_POUT;  /*  Duration  of  P  out  Transient  */ 
float  INIT_PRDEN;  /*  Initial  Value  of  Power  Density  */ 
float  FIN_PRDEN;  f*  Final  Value  of  Power  Density  */ 
float  DUR_PRDEN;  /*  Duration  of  Power  Density  Transient  */ 
float  t_DELAY;  /*  Time  Delay  prior  to  Transient  */ 
float  t_ AFTER;  f*  Window  of  Time  After  Transient  */ 

float  delta_t;  /*  Time  Step  */ 

float  INIT_W;  /*  Initial  Guess  at  Mass  Flowrate  */ 

); 

void  showdefaultsf  struct  fuelelem  *e_ptr ); 

void  fuelelemsum(  struct  fuelelem  *e_ptr,  struct  fuelelem  *f_ptr ); 

void  showfuelpart  (  struct  fuelpart  *g  ptr ); 

void  showconditions  (  struct  conditions  *h_ptr ); 

double  h2density(double  xx,  double  yy); 

double  h2viscosity(double  xx); 

doul  .e  h2heatxfer(doubIe  xx); 

double  h2cp(double  xx); 

float  h2prandle(float  xx); 

float  h2H(float  xx); 

I*  ***  BEGINNING  OF  COMPUTATIONAL  CODE  *******************  */ 
main() 

{ 

/*  Define  variables:  */ 

int  ch,  chr, 
float  Wl,  W2; 
float  W; 
double  Werror, 
float  pinl,  pinll,  pinin,  pout; 
float  pbarl,  pbarll,  pbarlll; 
float  uin,  uout,  ubarll; 

float  Tin,  TinI,  Tinll,  TinlH,  Tout,  Tbarl,  Tbarll,  TbarHI; 

float  Pinl,  Pinll,  Pinlll,  Pout,  Pbarl,  Pbarll,  Pbarlll; 

float  Pinll,  Pinlll,  Pinin  1,  Poutl,  Pbarll,  Pbarlll,  Pbarini; 

float  PDROPI1,  PDROPm,  PDROPIH1,  PDROPI2,  PDROPII2,  PDROPHI2; 

float  Hinl, Hinll,  HinUI,  Hinim,  Houtl; 

float  Hin2,  HinI2,  HinII2,  HinIII2,  Hout2,  Htemp; 

double  Herror, 

float  Deq_in,  Deq_out; 

float  AM,  AoutI,  ACF,  Ainll,  Aoutll,  AHF,  Ainlll,  Aoutlll; 

float  Abarll; 

float  voll,  volll,  volIII; 

float  RI,  RD,  Rin,  RTOTAL; 

float  Rif,  RIa,  RIk,  Rim,  Rlla; 

float  RCF,  Rel,  Rell,  ReUI; 

float  RUlf,  RHIa,  Rink,  Rlllm,  mfout; 

float  Qbed,  powerdensity,  vout,  power,  powerout; 

double  time,  transtime,  deltat; 

float  mal,  ma2,  ma3,  ma4,  mat; 

float  U 1,  U2,  U 3,  U 4,  UT,  UBAR,  Cpbar; 

float  fl,  f2,  f3.  fbar, 
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float  Nu,  Pr.  k,  h; 

float  II,  12, 13,  Ml,  M2,  M3; 

float  AA,  BB,  DD,  EE,  N; 

float  Vcv,  Av; 

float  Tf2,  Tfl,  Qs; 

float  WA,pfuelave; 

static  struct  fuelelem  A  = 

{  .04458,  .02936,  .00187,  .0000027,  .685,  .012432,  .0005, 

.4,  .00076,  .23,  .95, 1.1,  .265 

}; 

static  struct  fuelelem  B  = 

{  .04458,  .02936,  .00187,  .0000027,  .685,  .012432,  .0005, 

.4,  .00076,  .23,  .95, 1.1,  .265 

}; 

static  struct  fuelpart  C  = 

{  "UC2",  "Low  D  Car”,  "High  D  Car",  ”ZrC",  .0001 17,  .00015, \ 

.0002,  .00025, 10500, 1000,  1900,  6300, 200,  3000,  3000,  \ 

200,30, 1.5,3,40 

}; 

static  struct  conditions  D  = 

{  300, 300, 1, 2000, 2000, 1, 1915, 1915, 1, 0,  2,  1, 1,  8,  .01,  .1 

1; 

FILE  *out; 

out  =  fopen  (  "c:power.dat","w+"  ); 
ch  =  1; 

while  (ch  =  1 ) 

{ 

showdefaults  ( &A ); 

printf('ViCHANGE  ANOTHER  PARAMETER  FOR  ELEMENT  A?  (1=Y  or  0=N)"); 
scanf("%i",  &ch); 

} 

fuelelemsum  ( &A,  &B  ); 
ch  =  1; 

while  (ch  =  1 ) 

{ 

showfuelpart  (  &C ); 

printf('ViCHANGE  ANOTHER  PARAMETER?  (1=Y  or  0=N)”); 
scanf("%i",  &ch); 

} 

ch  =  1; 

while  (ch  =  1  ) 

{ 

showconditions  ( &D ); 

printf("\nCHANGE  ANOTHER  PARAMETER?  (1=Y  or  0=N)"); 
scanf(”%i",  &ch); 

1 

printf("Data  Sample  Frequency?(print  every  xth  data  point.)Nn"); 
scanf("%f ’,  &EE); 

fprintf(  out," VALUES  FOR  THE  FOLLOWING  INITIAL  CONDITIONS^"); 
fprintf(  out,"Vn  Initial  Inlet  Temp  =  %.2f  K\n",  D.INIT_TIN); 
fprintf(  out,"  Initial  Inlet  Press  =  %.2f  kPaNn",  D.INIT_PIN); 
fprintf(  out,"  Initial  Outlet  Press  =  %.2f  kPaViVi",  D.INIT_POUT); 
fprintf(out,  "timeNttr  t\tHI2StNtHII2NNtHIII2NMlout2\tpdenNtW2\n"); 
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W1  =  D.INITJW; 

AinI  =  PI  *  (A.OUTRAD  *  A.OUTRAD  -  A.CFRAD  *  A.CFRAD); 

Aoud  =  2  *  PI  *  A.CFRAD  *  A.LENGTH; 

ACF  =  Aoud  *  A.CFPOR; 

Ainn  =  2  *  PI  *  (A.CFRAD  -  A.CFTHK)  *  A.LENGTH  *  A.PBEDVOID; 

AoudI  =  2  *  PI  *  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK)  *  A.LENGTH  *  \ 

A.PBEDVOID; 

Abarll  =  .5  *  (Ainll  +  Aoutll); 

AinlD  =  2  *  PI  *  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK  -  A.HFTHK)  *  \ 

A.LENGTH; 

AoutUI  =  PI  *  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK  -  A.HFTHK)  *  \ 

(A.CFRAD  -  A.CFTHK  -  A.PBEDTHK  -  A.HFTHK); 

AHF  =  AinlH  *  A.HFPOR; 

Deq_in  =  4  *  AinI  /(2  *  PI  *  (A.OUTRAD  +  A.CFRAD)); 

Deq_out  =  4*  AoutHI  /  (2  *  PI  *  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK  -  \ 

A.HFTHK)); 

voU  =  AinI*  A.LENGTH+ACF*  A.CFTHK; 

volH  =  A.PBEDVOID*PI*((A.CFRAD-A.CFTHK)*(A.CFRAD-A.CFTHKHA.CFRAD-  \ 

A.CFTHK-A.PBEDTHK)*(A.CFRAD-A.CFTHK-A.PBEDTHK))*A.LENGTH; 
volffl  =  Aoutm*A.LENGTH  +  AHF*  A.HFTHK; 
transtime=0; 

N  =  0; 

for  (time  =  0;  time  <=  D.t_DELAY+D.delta_t;  time  +=  D.delta_t) 

l 

powerdensity  =  D.INIT_PRDEN; 

W  =  Wl; 

Qbed  =  powerdensity  *  1000000000  *  PI  *  ((A.CFRAD  -  A.CFTHK)  *  \ 

(A.CFRAD  -  A.CFTHK)  -  (A.CFRAD  -  A.CFTHK  -  A.PBEDTHK)  *  \ 
(A.CFRAD  -  A.CFTHK  -  A.PBEDTHK))  *  A.LENGTH; 

Hinll  =  h2H(D.INIT_TIN); 

PinI  =  D.INITPIN; 

Pinll  =  D.INIT_PIN; 

Pinm  =  DJNIT_POUT; 

Pout  =  D.INIT_POUT; 

Tin  =  D.INITTIN; 

Tout  =  Tin; 
do 
{ 

W- Wl; 

Houtl  =  Qbed/W  + Hinll; 
do 
{ 

Htemp  =  h2H(Tout); 

Herror  =  Houtl  -  Htemp; 

Tout  =  Tout  +  Herror/1 00000; 

}  while  (  Herror/ 10000  >  .001 ); 

/*  Set  up  all  variables  */ 

Tbarll  =  .5  *  (Tin  +  Tout); 

Pbarl  =  .5  *  (PinI  +  Pinll); 

Pbarll  =  .5  *  (Pinll  +  Pinhl); 

Pbarlll  =  .5  *  (PinlD  +  Pout); 
pbarl  =  h2density(Tin,  Pbarl); 
pbarll  =  h2density (Tbarll,  Pbarll); 
pbarIII=  h2density(Tout,  Pbarlll); 
uin  =  h2viscosity(Tin); 
ubarll  =  h2viscosity(TbarII); 
uout  =  h2viscosity(Tout); 
pinl  =  h2density(Tin,  PinI); 


126 


pinll  =  h2density(Tin,  Pinll); 
pinlll  =  h2density(Tout,  PinUI); 
pout  =  h2density(Tout,  Pout); 

/*  Inlet  Resistance  Calculations  */ 

Rel  =  W*Deq_in  /  (uin*AinI); 

Rif  =  ,138*pow(ReI,  15 1)*(  A.LENGTH/Deq_in)  /  (2*pbarI*AinI*AinI); 

RIa  =  (l/(Ainn*AinH)  -  l/(AinI*AinI))  /  (2*pbarl); 

RIk  =  ( 1  -( AinI/AoutI))*(  1  -(Ainl/AoutI))  /  (2*pbarI*AinI*AinI)  \ 

+  .5  *  ( 1  -( ACF/AoutI))*(  1  -( ACF/AoutI))  /  (2*pinII*ACF*ACF)  \ 

+  .5  *  (1  -( AMI/ACF))*(1 -( Ainll/ACF))  /  (2*pinn*AinD*AinII); 

Rim  =  A.MFIN  /  (pbarI*AinI*AinI); 

Rim  =  Rim  +  A.MFIN/(pbarI*AoutI*AoutI); 

RCF  =  150*uin*(l-A.CFPOR)*(l-A.CFPOR)/(A.CFDp*A.CFDp*pinII*AoutI*W*  \ 
A.CFPOR*  A.CFPOR*  A.CFPOR)  \ 

+  1.75*(l-A.CFPOR)/(A.CFDp*pinII*AoutI*AoutI*  A.CFPOR*  A.CFPOR”'  \ 
A.CFPOR); 

RCF  =  RCF  *  A.CFTHK; 

RI  =  RIf/2  +  RIa  +  RIk  +  RIm  +  RCF; 

/*  Resistance  Calculations  for  the  Particle  Bed  */ 

RH  =  150*ubarn*(l-A.PBEDVOID)*(l-A.PBEDVOID)/(A.PBEDDp*A.PBEDDp*\ 
pbaxn*AbarH*  W*  A.PBEDVOID*  A.PBEDVOID)  \ 

+  1.75*(l-A.PBEDVOID)/(A.PBEDDp*pbarn*Abarn*Abarn*A.PBEDVOID); 
RII  =  RII  *  A.PBEDTHK; 

RHa  =  ( l/(pinni*Ainin)  -  l/(pinH*Ainn)  )*((AMI+Ainm)/  \ 

(2*Ainn*Ainm)); 

RII  =  RH  +  RHa; 

/*  Resistance  for  Exit  Plenum  */ 

Relll  =  W*Deq_out  /  (uout*Aoutm); 

RHIf  =  .138*pow(Rem,  -.151)*(A.LENGTH/Deq_out)/(2*pbarin*Aoutm\ 

*AoutIII); 

RHIa  =  (l/(Aoutm*Aoutm>l/(Aoutn*Aoutn»/(2*pbarin); 

Rlllk  =  .5*(l-AHF/Aoutn)*(l-AHF/Aoutn)/(2*pinm*AHF*AHF)\ 

+  (l-AHF/Ainm)*(l-AHF/Ainm)/(2*pihffl*AHF*AHF)  \ 

+  (l-Aoutm/Ainffl)*(l-Aoutni/Ainm)/(2*pbaini*  \ 

Aoutin*Aoutin); 
mfout  =  A.MFOUT; 

Rnim  =  mfout  /  (pbarin*Aoutin*Aoutin)  +  mfout  /  (pbarlll  \ 

*Ainin*Ainin); 

Rin  =  Rfflf/2  +  RHIa  +  Rlllk  +  Rnim; 

RTOTAL  =  RI  +  RH  +  RIH; 

W1  =  (PinI  -  Pout)*  lOOO/RTOTAL; 

W1  =sqrt(Wl); 

Werror  =  W  -  Wl; 

Werror  =  fabs(Werror); 

PDROPI1=RI*W1*W1/1000; 

PDROPII1=RII*W1*W1/1000; 

PDROPim=Rni*Wl*Wl/1000; 

Pinn=PinI-PDROPIl; 

Pinm=Pout+PDROPffll; 

PbarI=.5*(PinI+Pinn); 

PbarII=.5*(PinII+Pinni); 

Pbarin=.5*(Pinin+Pout); 
vout  =  Wl/(pout*Aoutin); 
power  =  powerdensity; 

)  whilefWerror  >  .0000001); 
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HinI2=HinIl; 

HinIIl=HinIl; 

HinII2=HinIIl; 

Hinnil=Houtl; 

Hinni2=Hinini; 

Hout2=Houtl; 

powerout  =  W*(Hout2-HinI2); 

DD  =  finod(N,  EE); 
if  (DD  ==  0) 

{ 

W1=W1*1000; 

fprintf(out,''%Af\t%N%fSt%f\t%f\t%f4%f't%f^%f\t%N%f\t%<N«'',\ 
time,  power,  Wl,  Tin,  Tbarll,  Tout,  powerout,  \ 

PinI,  Pinll,  PinUI,  Pout,  vout); 

printf("time=%f''tTin=%fStTout=%NW=%f\n",  time,  Tin,  Tout,  Wl); 
W1=W1/1000; 

} 

DD  =  0; 

N  +=  1.0; 
if  (Tout  >  2700) 

goto  meltdown; 


} 

Pinll  =  PinI; 

Pinll  1  =  Pinll; 

PinlTll  =  Pinlil; 

Poutl  =  Pout; 

Pbarll  =  Pbarl; 

Pbarlll  =  Pbarll; 

Pbarllll  =  Pbarlll; 

mal  =  ((C.L3RAD*C.L3RAD*C.L3RAD)-(C.L2RAD*C.L2RAD*C.L2RAD))*C.L3DEN/\ 
(3*C.L3RAD*C.L3RAD); 

ma2  =  ((C.L2RAD*C.L2RAD*C.L2RADHC.LlRAD*C.LlRAD*C.LlRAD))*C.L2DEN/\ 
(3*C.L3RAD*C.L3RAD); 
ma3  =  ((C.L1RAD*C.L1RAD*C.L1RAD)- 
(C.FUELRAD*C.FUELRAD*C.FUELRAD))*C.L  1DEN/  \ 

(3  *C.L3RAD*C.L3RAD); 

ma4  =  (C.FUELRAD*C.FUELRAD*C.FXJELRAD)*C.FUELDEN  /  (3*C.L3RAD*C.L3RAD); 
mat  =  mal+ma2+ma3+ma4; 

pfuelave=(l/(  1.3333*PI*pow(C.L3RAD,3)  )  )*(  C.FUELDEN*pow(C.FUELRAD3)\ 

+  C.L1DEN*(  pow(C.LlRAD,3)-pow(C.FUELRAD,3) )  +  C.L2DEN*(  pow(C.L2RAD,3)\ 
-pow(C.LlRAD,3) )  +  C.L3DEN*(  pow(C.L3RAD,3)-pow(C.L2RAD,3) ) ); 

Cpbar  =  (mal*C.L3Cp  +  ma2*C.L2Cp  +  ma3*C.LlCp  +  ma4*C.f?UELCp)^iat; 

U1  =  C.L3RAD*C.L2RAD*C.L3K/(C.L3RAD*C.L3RAD*(C.L3RAD-C.L2RAD)); 

U2  =  C.L2RAD*C.L1RAD*C.L2K/(C.L3RAD*C.L3RAD*(C.L2RAD-C.L1RAD)); 

U3  =  C.L1RAD*C.FUELRAD*C.L1K/(C.L3RAD*C.L3RAD*(C.L1RAD-C.FUELRAD)); 

U4  =  2*C.FUELK/(3*C.L3RAD); 

UT  =  1/C1/U1  +  1/U2  +  1/U3  +  1/U4); 
fl  =  UT/U1; 
f2  =  UT/U2  +  fl; 

G  =  UT/U3  +  £2; 

fbar  =  (mal*C.L3Cp*fl  +  ma2*C.L2Cp*(fl+G)  +  ma3*C.LlCp*(f2+G)  +\ 
ma4*C.FUELCp*(G+l))/(2*mat*Cpbar); 

V  cv=volII/A  .PB  ED  V  OID ; 

Av=6*(l-A.PBEDVOID)/A>PBEDDp; 

TinI=Tin; 

Tfl  =  Tbarll; 

TinII=TinI; 
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TinIII=Tout; 

deltat=D.delta_t; 

N  =  0.0; 

for  (time=D.t_DELAY ;  time  <  (D.t_DELAY+D.DUR_PRDEN+D.t_ AFTER);  \ 
time  +=  deltat) 

{ 

poweidensity  =  D.INIT_PRDEN  +  ((D.FIN_PRDEN  -  D.INIT_PRDEN)  \ 

*  transtime  /  D.DUR_PRDEN); 
if  (poweidensity  >=  D.FIN_PRDEN) 
powerdensity  =  D.FIN_PRDEN; 
k  =  h2heatxfer(TbarII); 

Pr  =  h2prandle(TbarII); 

ReH  =  (W*6)/(AbarI]*A.PBEDVOID*ubarn*Av); 

Nu  =  .8*pow(ReII,  ,7)*pow(Pr,  .33); 
h  =  Nu*k*2/A.PBEDTHK; 

UBAR  =  UT*ta/(fbar*h+UT); 
power  =  powerdensity; 

Qs  =  power*  1000000000/A  v; 

BB=mat*Cpbar*Vcv*Av; 

AA=UBAR*Vcv*Av; 
if  (transtime  —  0); 

Tfl  =Qs*  V  cv*  Av/(AA)  +  Tbarll; 
if  (transtime  >  1) 

{ 

if  (D.DUR_TIN  >  0) 

{ 

TinI=D.INIT_TIN+(D.FIN_TIN  -  D.INIT_TIN)*(transtime/D.DlJR_TIN); 
HinI2=h2H(TinI); 

} 

if  (transtime  >=  D.DUR_TIN) 

{ 

TinI  =  D.FINJITN; 

HinI2  =  h2H(TinI); 

} 

} 

if  (D.DUR_PIN  >  0) 

Pinl  1  =D.INTT_PIN+(D.FIN_PIN-D.INIT_PIN)*(transtime/D.DlJR_PIN); 
if  (transtime  >=  D.DUR_PIN) 

Pinll  =  D.FIN_PIN; 
if  (D.DUR_POUT  >  0) 

Poutl=D.INIT_POUT+(D.FrN_POUT-D.INrr_POlJT)*(transtime/D.DlJR_POUT); 
if  (transtime  >=  D.DUR_POUT) 

Poutl  =  D.FIN_POUT; 

/*  Control  Volume  I  */ 

Ml  =pinI*volI; 

11  =  ,5*A.LENGTH/AinI+  (A.OUTRAD-A.CFRAD)/AoutI  +  A.CFTHK/ACF, 

W2=(W*I  1  +deltat*(PinI  1-PinII  1  )*  1 000)/(1 1 4deltat*RI*  W); 

Hinl32  =  (HinHl  +  (deltat/Ml)*((von*(PbarIl-PbarI)*1000/deltat)  +\ 

W*HinI2  +  0.00*Qbed))/(l+deltat*W/Ml); 

/*  Control  Volume  n  */ 

Tf2  =  (Tfl  +  deltat*(Qs*Vcv*Av  +  AA*Tbarn)/BB)/(  1  +  deltat*  AA/BB); 

Qbed  =  AA*(Tf2-Tbarn); 

Tfl  =Tf2; 

M2  =  (pbar!I*voin  +  pfuelave*voin*(l-A.PBEDVOID))/2; 

12  =  A-PBEDTHK/AbarH; 

W2=(  W*I2+deltat*(Pinn  l-PinlH 1  )*  1 000)/(I2+deltat*RII*  W); 

Hinm2  =  (Hinmi  +(deltat/M2)*((Vcv*(Pbarni-Pbarn)*1000/deltat)\ 

+  W*HinII2  +  1.00*Qbed  »/(l+deltat*W/M2); 
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/*  Control  Volume  in  */ 

M3  =  pinIH*volffl; 

13  =  A.HFTHK/AoutII  +  (A.CFRAD-A.CFTHK-A.PBEDTHK  \ 
-A.HFTHK)/Ainm  +  ,5*A.LENGTH/AoutIII; 

W2=(  W*I3+deltat*(Pinini  -Pout 1 )*  1 000)/(I3+deltat*RIII*  W); 

Hout2  =  (Houtl  +  (deltat/M3)*(  (volIII*(Pbarnil-Pbarni)*1000/deltat)\ 
+  W*HinIII2  +  0.00*Qbed  ))/(l+deltat*W/M3); 


I*  ******************************************************************  * j 

Pbarl  =  Pbarll; 

Pbarll  =  Pbarlll; 

Pbarlll  =  Pbarini; 

PinI  =  Pinll; 

Pout  =  Poutl; 

PinU  =  Pinll  1; 

PinlTI  =  PinlTIl; 
do 
{ 

Htemp  =  h2H(TinI); 

Herror  =  HinI2  -  Htemp; 

TinI  =  TinI  +  Henor/100000; 

}  while  (  Henor/100000  >  .001 ); 

do 

{ 

Htemp  =  h2H(Tout); 

Heiror  =  Hout2  -  Htemp; 

Tout  =  Tout  +  Henor/100000; 

}  while  (  Henor/100000  >  .001  ); 
do 
{ 

Htemp  =  h2H(Tmfl); 

Herror  =  HinII2  -  Htemp; 

Tinll  =  Tinll  +  Henor/100000; 

}  while  ( Henor/100000  >  .001 ); 
do 
{ 

Htemp  =  h2H(TinIII); 

Herror  =  HinIII2  -  Htemp; 

Tinin  =  TinlH  +  Henor/100000; 

}  while  ( Herror/100000  >  .001 ); 


TbarI=.5*(TinI+TinII); 

TbarII=.5*(TinII+TinIII); 

TbarIII=.5*(TinIII+Tout); 

Hi  nil  =  HinI2; 

Hi  nil  1  =  HinII2; 

Hinllll  =  HinIII2; 

Houtl  =  Hout2; 
pinl  =  h2density(TinI,  PinI); 
pinll  =  h2density(TinII,  Pinll); 
pinlll  =  h2density(TinIII,  Pinlll); 
pout  =  h2density(Tout,  Pout); 
pbarl  =  h2density(TbarI,PbarI); 
pbarll  =  h2density(TbarII,  Pbarll); 
pbarll  =  h2density(Tout,  Pbarlll); 


uin  =  h2viscosity(TbarI); 
ubarll  =  h2viscosity(TbarII); 
uout  =  h2viscosity(TbarlII); 

/*  Inlet  Resistance  Calculations  */ 

Rel  =  (W*Deq_in)/(uin*  AinI); 

Rif  =  ,138*pow(ReI,  15 l)*(A.LENGTH/Deq_in)/(2*pbarI* AinI* AinI); 

RIa  =  (l/(AinH*AinII)  -  l/(AinI*AinI))/(2*pbarI); 

RIk  =  ( 1  -(AinI/AoutI))*(  1  -(Ainl/AoutI))  /  (2*pbari*  AinI*  AinI)  \ 

+  .5  *  (l-(ACF/AoutI))*(l-(ACF/AoutI))  /  (2*pinn*ACF*ACF)  \ 

+  .5  *  ( 1  -( Ainll/ ACF)) * ( 1  -( Ainll/ACF))  /  (2*pinn*Ainn*AinII); 

Rim  =  A.MFIN  /  (pbarl*  AinI*  AinI); 

Rim  =  Rim  +  A.MFIN/(pbarI*AoutI*AoutI); 

RCF  =  150*uin*(l-A.CFPOR)*(l-A.CFPOR)/(A.CFDp*A.CFDp*pinII*AoutI*W*  \ 
A.CFPOR*  A.CFPOR*  A.CFPOR)  \ 

+  1.75*(l-A.CFPOR)/(A.CFDp*pinII*AoutI*AoutI*  A.CFPOR*  A.CFPOR*  \ 
A.CFPOR); 

RCF  =  RCF  *  A.CFTHK; 

RI  =  RIf/2  +  RIa  +  RIk  +  Rim  +  RCF; 

/*  Resistance  Calculations  for  the  Particle  Bed  */ 

Rn  =  150*ubarH*(l-A.PBEDVOID)*(l-A.PBEDVOID)/(A.PBEDDp*A.PBEDDp*\ 
pbarn*AbarII*W*  A.PBEDVOID*  A.PBEDVOID)  \ 

+  1.75*(l-A.PBEDVOID)/(A.PBEDDp*pbarn*AbarII*Abarn*  A.PBEDVOID); 
Rn  =  Rn  *  A.PBEDTHK; 

Rlla  =  ( l/(pinm*Ainin)  -  l/(pinn*AMI)  )*((Ainn+Ainin)/  \ 

(2*Ainn*Ainin)); 

RH  =  RH  +  RHa; 

/*  Resistance  for  Exit  Plenum  */ 

Rein  =  W*Deq_out  /  (uout*Aoutni); 

Rim  =  .138*pow(Reffl,  -,151)*(A.LENGTH/Deq_out)/(2*pbarm*Aoutm\ 

*Aoutm); 

Rnia  =  (l/(Aoutin*Aoutm)-l/(Aoutn*Aoutn))/(2*pbarm); 

Rfflk  =  .5*(l-AHF/Aoutn)*(l-AHF/Aoutn)/(2*pinm*AHF*AHF)  \ 

+  (l-AHF/Ainm)*(l-AHF/Ainm)/(2*pinm*AHF*AHF)  \ 

+  (l-Aoutm/Ainm)*(l-Aoutm/Ainffl)/(2*pbarm*  \ 

AoutIII*Aoutm); 
mfout  =  A.MFOUT; 

RHIm  =  mfout  /  (pbarm*Aoutm*  AoutHI)  +  mfout  /  (pbarlll  \ 

*Ainm*Ainm); 

RIII  =  RHIf/2  +  RHIa  +  Rlllk  +  RHIm; 

/*  TOTALS  */ 

RTOTAL  =  RI  +  RII  +  RIU; 

W1  =  (PinI  -  Pout)*  lOOO/RTOTAL; 

Wl=sqrt(Wl); 

W  =  Wl; 

PDROPI1  =  RI*W*W/1000; 

PDROPII1  =  RII*W*W/1000; 

PDROPIH1  =Rm*W*W/1000; 

Pinlll  =  Pinll  -  PDROPI1; 

Pinffll  =  Poutl  +  PDROPmi; 

Pbarll  =  .5*(PinIl+Pinni); 

Pbarlll  =  ,5*(PinIIl+Pinmi); 

Pbarllll  =  ,5*(Pinmi+Poutl); 
vout  =  W/(pout*Aoutm); 
powerout  =  W*(Hout2-HinI2); 

DD  =  finod(N,  EE); 
if  (DD  ==  0) 

{ 
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W  =  W*1000; 

fprintf(out,"%.4f\t  %.4f\t  %.3f*  %.3f\t  %3H  %.2N.  %.2f\t  %.2H\ 

%.2H  %.2H  %.2N  %.3fV,  \ 

time,  power,  W,  TinI,  Tbarll,  Tout,  powerout,  \ 
Pinll,  PinUl,  Pinini,  Poutl,  vout,  Tfl); 
printf("time=%f\tTinI=%iNtTout=%NW=%fV,  time,  TinI,  Tout,  W); 

W  =  W/1000; 

} 

DD  =0; 

N  +=  1.0; 
if  (Tout  >2700) 

goto  meltdown; 
transtime  +=  deltat; 

1 

meltdown: 

if  (Tout  >2700) 

{  printfCTout  >  2700”); 

fjprintf(out,'Temp  exceeds  2700"); 

} 


/*  ******************************************************************  *j 

void  showdefaults(  stnict  fuelelem  *e_ptr ) 

{ 

int  num; 
float  entry; 
int  c; 

_clearscreen(0); 

printf('THE  FOLLOWING  ARE  PARAMETERS  FOR  PPBR  FUEL  ELEMENT :\n\n"); 
printf("l.  Outer  Radius  =  %f  mV,  e_ptr->OUTRAD); 

printf("2.  Cold  Frit  Radius  =  %f  mV,  e_ptr->CFRAD); 

printf(”3.  Cold  Frit  Thickness  =  %f  mV,  e_ptr->CFTHK); 
printf("4.  Dia  of  Cold  Frit  Particle  =  %.8f  m'vi",  e _ptr->CFDp); 
printf("5.  Cold  Frit  Porosity  =  %fV,  e_ptr->CFPOR); 
printf("6.  Particle  Bed  Thickness  =  %f  m\n",  e_ptr->PBEDTHK); 
printf("7.  Fuel  Particle  Diameter  =  %f  mV,  e_ptr->PBEDDp); 

printf("8.  Particle  Bed  Void  Fraction  =  %fV,  e_ptr->PBEDVOID); 
printf(”9.  Hot  Frit  Thickness  =  %f  mV,  e_ptr->HFTHK); 
printf("10.  Hot  Frit  Porosity  =  %fV,  e _ptr->HFPOR); 
printf('T  1.  Inlet  Manifold  Factor  =  %1V,  e_ptr->MFIN); 
printf("  12.  Exit  Manifold  Factor  =  %fV,  e_ptr->MFOUT); 

printfC  13.  Fuel  Element  Length  =  %f  mV,  e_ptr->LENGTH); 
printfCVnNnCHANGE  A  PARAMETER?  (1=Y  or  0=N)"); 
scanf(  "%i",  &c ); 
if  ( c  =  0  ) 
return; 

printfCVLEASE  ENTER  THE  PARAMETER  NUMBER  AND  VALUE  (eg  12, 1.705):V); 
scanf(  "%i,  %f ',  &num,  &entry  ); 
switch  (num) 

{ 

case  1: 

e_ptr->OUTRAD  =  entry; 
break; 
case  2: 

e_ptr->CFRAD  =  entry; 
break; 

case  3:  _ 

e_ptr->CFTHK  =  entry; 
break; 
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case  4: 


e_ptr->CFDp  =  entry; 
break; 

case  5: 

e_ptr->CFPOR  =  entry; 
break; 

case  6: 

e_ptr->PBEDTHK  =  entry; 
break; 

case  7: 

e_ptr->PBEDDp  =  entry; 
break; 
case  8: 

e_ptr->PBED  V  OID  =  entry; 
break; 

case  9: 

e_ptr->HFTHK  =  entry; 
break; 
case  10: 

e_ptr->HFPOR  =  entry; 
break; 
case  11: 

e_ptr->MFIN  =  entry; 
break; 
case  12: 

e_ptr->MFOUT  =  entry; 
break; 
case  13: 

e_ptr->LENGTH  =  entry; 
break; 
default: 

printf("\n\nTRY  ANOTHER  PARAMETERVi”); 

_clearscreen(0); 

break; 

} 

return; 

1 

jAt  ****+*************************************************************  Atf 

void  fuelelemsum  ( struct  fuelelem  *e_ptr,  struct  fuelelem  *f_ptr ) 

f 

clears  crcen  (0)* 

printf  ("PARAMETERW'tELEMENT  ANtNtBASEUNEnW); 
printf("l.  Outer  RadiusMNlNt%f\t\t%f  mNn",  e_ptr->OUTRAD,  f_ptr->OUTRAD); 
printf("2.  Cold  Frit  RadiusNt\t%fNt\t%f  mVi",  e_ptr->CFRAD,  f_ptr->CFRAD); 
printf("3.  Cold  Frit  ThkknessMWcNWcf  mV’,  e_ptr->CFTHK,  f_ptr->CFTHK); 
printf("4.  Dia  of  Cold  Frit  Particle\t%.8f\t\t%.8f  mNn”,  e_ptr->CFDp,  f_ptr->CFDp); 
printf("5.  Cold  Frit  Porosity\t\t%fNtM%f\n",  e_ptr->CFPOR,  f_ptr->CFPOR); 
printf("6.  Particle  Bed  thickness\t%f^t%f  mNn",  e_ptr->PBEDTHK,  f_ptr->PBEDTHK); 
printf("7.  Fuel  Particle  DiameterNt%f\fNt%f  mNn”,  e  j>tr->PBEDDp,  f_ptr->PBEDDp); 
printf("8.  Particle  Bed  Void  FractionNt%fNtNt%fNn”,  e_ptr->PBEDVOID,  f_ptr->PBEDVOID); 
printf("9.  Hot  Frit  ThicknessNt\t%fNtM%f  mNn",  e_ptr->HFTHK,  f_ptr->HFTHK); 
printf("10.Hot  Frit  PorosityNfst%fNtNt%fn”,  e_ptr->HFPOR,  f_ptr->HFPOR); 
print f(”  11. Inlet  Manifold  Factorit%fNtNt%fNn",  e_ptr->MFIN,  f_ptr->MFIN); 
printf("12.Exit  Manifold  FactoiNt\t%fNtNt%fn”,  e_ptr->MFOUT,  f_ptr->MFOUT); 
printf("13.Length  of  Element\tNt%f  mMNt%f\n",  e_ptr->LENGTH,  f_ptr->LENGTH); 
printf("NriNnHIT  ANY  KEY  TO  CONTINUE"); 
getchQ; 


y*  ******************************************************************  *y 

void  showfuelpait  {  struct  fuelpart  *g_ptr ) 

{ 

int  num; 
float  entry; 
intc; 

clearscieen(O); 

printf('THE  FOLLOWING  ARE  PARAMETERS  FOR  THE  FUEL  PARTICLESn"); 
printf("l.  Fuel  Material  =  %sNn",  g_ptr->FUELMAT); 
printf("2.  Layer  1  Material  =  %s\n",  g_ptr->LlMAT); 
printf(”3.  Layer  2  Material  =  %sNn",  g_ptr->L2MAT); 
printf("4.  Layer  3  Material  =  %sNn",  g_ptr->L3MAT); 
printf("5.  Radius  of  Fuel  =  %f  m\n",  g_ptr->FUELRAD); 
printf("6.  Radius  of  Layer  1  =  %{  mV,  g_ptr->LlRAD); 
printf("7.  Radius  of  Layer  2  =  %f  m\n",  g_ptr->L2RAD); 
printf("8.  Radius  of  Layer  3  =  %f  m\n",  g_ptr->L3RAD); 
printf("9.  Density  of  Fuel  =  %f  kg/m3\n",  g_ptr->FUELDEN); 
printf("10.Density  of  Layer  1  =  %f  kg/m3\n",  g_ptr->LlDEN); 
printf("  11. Density  of  Layer  2  =  %fkg/m3Nn",  g_ptr->L2DEN); 
printf("12.Density  of  Layer  3  =  %f  kg/m3\n",  g_ptr->L3DEN); 
printf("13.  Cp  of  Fuel  =  %f  J/kgKNn",  g_ptr->FUELCp); 

printf("14.  Cp  of  Layer  1  =  %f  J^cgKSn",  g_ptr->LlCp); 
printf("15.  Cp  of  Layer  2  =  %f  J/kgKNn",  g_ptr->L2Cp); 
printf("16.  Cp  of  Layer  3  =  %f  J/kgKNn",  g_ptr->L3Cp); 
printf('T7.  k  for  Fuel  =  %f  W/m2KNn",  g_ptr->FUELK); 
printf("18.  k  for  Layer  1  =  %f  W/m2KNn",  g_ptr->LlK); 
printf("  19.  k  for  Layer  2  =  %f  W/m2KNn”,  g_ptr->L2K); 
printf("20.  k  for  Layer  3  =  %f  W/m2K>o",  g  ptr->L3KV. 
printf("NnCHANGE  A  PARAMETER?  (1=Y  or  0=N)”); 
scanf(  "%i",  &c ); 
if  ( c  ==  0 ) 
return; 

printf(’NnENTER  PARAMETER  NUMBER  (  NOT  1-4  ),  VALUE  (eg  20, 40)W'); 
scanf(  "%i,  %f &num,  &entry  ); 
switch  (num) 

{ 

case  5; 

g_ptr-  >FUELR  AD  =  entry; 
break; 

case  6; 

g_ptr->LlRAD  =  entry; 
break; 

case  7: 

g_ptr->L2RAD  =  entry; 
break; 
case  8: 

g_ptr->L3RAD  =  entry; 
break; 

case  9: 

g_ptr->FUELDEN  =  entry; 
break; 
case  10: 

g_ptr->LlDEN  =  entry; 
break; 
case  11; 

g_ptr->L2DEN  =  entry; 
break; 
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case  12: 

g_ptr->L3DEN  =  entry; 
break; 
case  13: 

g_ptr->FUELCp  =  entry; 
break; 
case  14: 

g_ptr->LlCp  =  entry; 
break; 
case  15: 

g_ptr->L2Cp  =  entry; 
break; 
case  16: 

g_ptr->L3Cp  =  entry; 
break; 
case  17: 

e  ptr->FUELK  =  entry; 
break; 
case  18: 

g  ptr->L!K  =  entry; 
break; 
case  19: 

g_ptr->L2K  =  entry; 
break; 
case  20: 

g  ptr->L3K  =  entry; 
break; 
default: 

printf("CANNOT  CHANGE  THIS  PARAMETER...SORRY.V); 

getch(); 

break; 

} 

return; 

I 

/*  ******************************************************************  * j 

void  showconditions  (  struct  conditions  *h_ptr ) 

( 

int  num; 
float  entry; 
int  c; 

_clearscreen(0); 

printfCTHE  FOLLOWING  INITIAL  AND  FINAL  CONDITIONS  ARE  SETyiNn"); 
printf("  1 .  Initial  Inlet  Temperature  =  %f  K\n",  h_ptr->INIT_TIN); 

printf("2.  Final  Inlet  Temperature  =  %f  KNn",  h_ptr->FIN_TIN); 
printf("3.  Duration  of  Temperature  Change  =  %f  secNn",  h_ptr->DUR_TIN); 
printf("4.  Initial  Inlet  Pressure  =  %f  KPaW,  h_ptr->INIT_PIN); 

printf("5.  Final  Inlet  Pressure  =  %f  KPaNn”,  h_ptr->FIN_PIN); 

printf("6.  Duration  of  Pressure  Change  =  %f  secNn",  h_ptr->DUR_PIN); 
printf(”7.  Initial  Outlet  Pressure  =  %f  KPaW,  h_j)tr->INIT_POUT); 

printf("8.  Final  Outlet  Pressure  =  %f  KPaNo",  h_ptr->FIN_POUT); 

printf(”9.  Duration  of  Pressure  Change  =  %f  secNn",  h_ptr->DUR_POUT); 
printf('T0.  Initial  Power  Density  =  %f  GW/m3 Nn",  h_ptr->INIT_PRDEN); 

printf(”  1 1 .  Final  Power  Density  =  %f  GW/m3\n'\  h_ptr->FIN_PRDEN); 

printf("12.Duration  of  Power  Density  Change  =  %f  secNn",  h_ptr->DUR_PRDEN); 
printf("13.  Time  Delay  to  Transient  =  %f  secNn",  h_ptr->t_DELAY); 
printf("14.  Time  After  Transient  =  %f  secNn",  h_ptr->t_AFTER); 
printf("15.  Time  Step  =  %f  secNn",  h_ptr->delta_t); 

printf("16.  Initial  Guess  at  Mass  Flowrate  =  %f  kg/secNn",  h_ptr->INIT_W); 
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printf("NnCHANGE  A  PARAMETER?  (]  =  Y  or  0=N)"); 
scanf(  "%i",  &c ); 
if  ( c  =  0  ) 
return; 

printf(”ViENTER  PARAMETER  NUMBER,  VALUE  (eg  20,  40)\n"); 
scanf(  "%i,  %f",  &num,  &entry  ); 
switch  (num) 

{ 

case  1: 

h_ptr  >INTT_TIN  =  entry; 
break; 

case  2: 

h_ptr->HN_TIN  =  entry; 
break; 

case  3: 

h_ptr->DUR_TIN  =  entry; 
break; 

case  4: 

h_ptr->INIT_PIN  =  entry; 
break; 

case  5; 

h_ptr->HN_PIN  =  entry; 
break; 

case  6; 

h_ptr->DUR_PIN  =  entry; 
break; 

case  7: 

h_ptr->INTT_POUT  =  entry; 
break; 

case  8; 

h_ptr->FIN_POUT  =  entry; 
break; 

case  9; 

h_ptr->DUR_POUT  =  entry; 
break; 
case  10: 

h_ptr->INIT_PRDEN  =  entry; 
break; 
case  11: 

h_ptr->FIN_PRDEN  =  entry; 
break; 
case  12: 

h_ptr->DUR_PRDEN  =  entry; 
break; 

case  l.V 

h_ptr->t_DELAY  =  entry; 
break; 
case  14: 

h_ptr->t_AFTER  =  entry; 
break; 
case  15: 

h_ptr->delta_t  =  entry; 
break; 
case  16: 

h_ptr->INIT_W  =  entry; 
break; 
default: 

printfO'CANNOT  CHANGE  THIS  PARAMETER.. .SORRY.V); 
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getch(); 

break; 

} 

ret'im; 

} 

I *  ******************************************************************  *j 

double  h2density(double  xx,  double  yy) 

{ 

float  dens; 

dens  =  .08185  *  300  /  xx  *  yy  /  101 .325; 
return  dens; 

} 

! *  ******************************************************************  *i 

double  h2viscosity(double  xx) 

{ 

float  vise; 

vise  =  .000008963  *  pow(  xx/300,  .6733); 
return  vise; 

} 

I*  ******************************************************************  * j 

double  h2heatxfer(double  xx) 

{ 

float  teml; 
int  i; 

float  tem2; 

float  heat[45]  =  {0..0362,  .0665,  .0981,  .1282,  .1561,  .182,  .206,  ,228,\ 

.251,  .272,  .292,  .315,  .333,  .351,  .3665,  .384,  .398, \ 

.412,  .426,  .44,  .452,  .464,  .476,  .488,  .500,  .512,  \ 

.524,  .536,  .548,  .560,  .572,  .584,  .596,  .608,  .62,  \ 

.632,  .644,  .656,  .668,  .680,  .692,  .704,  .716,  .728  ); 

tem2  =  xx  /  50.0; 
i  =  floor(tem2); 

teml  =  heat[i]  +  (heat[i+l]  -  heat[i])  *  ((xx  -  (i  *  50))  /  50 ); 
return  teml; 

} 

******************************************************************  *^ 
double  h2cp(double  xx) 

{ 

float  temcp; 
if  ( xx  <=  420  ) 

temcp  =  4.1868  *  1000  *  (1.5395  +  .0150825  *  xx  -  4.02449e-5  *\ 
xx  *  xx  +  3.63544e-8  *  xx  *  xx  *  xx); 

else 

temcp  =  4.1868  *  1000  *  (3.58927  -  5.55096e-4  *  xx  +  6.94235e-7  *  xx\ 
*  xx  -  1 .45 1 55e-10  *  xx  *  xx  *  xx); 

return  temcp; 

} 

******************************************************************  tfj 

float  h2H(float  xx) 

{ 

float  temH; 
float  temHl; 
temH=0; 
temH  1=0; 
if  (xx  <=  420) 

{ 

temH  =4. 1868*1 000*(  1 .5395*xx  +  .0150825*xx*xx/2  -  4.02449e-5*  \ 

xx*xx*xx/3  +  3.63544e-8*xx*xx*xx*xx/4); 
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temH=temH  -  4.1868*1000*(1.5395*:>0  +  .0150825*50*50/2  -  4.02449e-5*  \ 
50*50*50/3  +  3.63544e-8*50*50*50*50/4); 

) 

else 

{ 

temH  1 =4 . 1 868*  1 000*(3 .58927 *xx  -  5.55096e-4*xx*xx/2  +  6.94235e-7*  \ 

xx*xx*xx/3  -  1.45155e-10*xx*xx*xx*xx/4); 
temHl=temHl  -  4.1868*1000*(3.58927*420  -  5.55096e-4*420*420/2  +  \ 

6.94235e-7*420*420*420/3  -  1.45155e-10*420*420*420*420/4); 
temH=temHl  +  4.1868*1000*(1.5395*420+.0150825*420*420/2-4.02449e-5*  \ 
420*420*420/3  +  3.63544e-8*420*420*420*420/4); 
temH=temH  -  4.1868*1000*(1.5395*50  +  .0150825*50*50/2  -  4.02449e-5*  \ 
50*50*50/3  +  3.63544e-8*50*50*50*50/4); 

} 

return  temH; 

} 

j*  ******************************************************************  * i 

float  h2prandle(float  xx) 

{ 

float  teml; 
int  i; 

float  tem2; 

float  pran[23]  =  { .000,  .712,  .719,  .706,  .690,  .675,  .664,  \ 

.659,  .664,  .676,  .686,  .703,  .715,  .733,  \ 

.748,  .763,  .778,  .793,  .808,  .823,  .838,  \ 

.853,  .868}; 

tem2  =  xx  /  100.0; 
i  =  floor(tem2); 

teml  =  pranfi]  +  (pran[i+l]  -  pranfi])  *  ((xx  -  (i  *  100))  / 100 ); 
return  teml; 

} 

/*  ******************************************************************  *  j 
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APPENDIX  E: 


TABLE  OF  SYMBOLS  AND  ACRONYMS 

A  cross  sectional  area  perpendicular  to  flow 
Av  particle  surface  area  per  unit  volume 

A^ii  smaller  areas  involved  in  expansion/contraction 
Cp  average  particle  specific  heat 
dp  effective  particle  diameter 
DP  particle  diameter 
D  equivalent  channel  diameter 

e  void  fraction  or  porosity 
Fa  equivalent  pressure  drop  due  to  acceleration 

(Fa)pBED  spatial  acceleration  losses  in  the  particle  bed 

Feq  equivalent  resistive  pressure  drop 
Ff  equivalent  pressure  drop  due  to  friction 
Fm  equivalent  pressure  drop  due  to  the  manifold  effect 
Fpbed  equivalent  pressure  drop  due  particle  bed  effects 
Fk  equivalent  loss  due  to  sudden  expansion/contraction 
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fj  friction  factor 

h  heat  transfer  coefficient 
H  specific  enthalpy 
lev  control  volume  inertia 

k  coolant  heat  conductivity 
KkJCeyKc  expansion/contraction  coefficients 

Lb  average  length  of  travel  through  particle  bed 
Lp  average  length  of  travel  through  plenums 
ma  particle  mass  per  unit  surface  area 

M  mass  within  the  control  volume 
Nup  particle  bed  Nusselt  number 

wetted  perimeter 

P  average  pressure  within  the  control  volume 
A Pay  pressure  drop  across  control  volume  due  to  the  resistance 

q  heat  energy 

q  ”  heat  generation  per  surface  area 
Q,  interphase  heat  transfer 

Re  Reynolds  Number 
RJo,at  total  equivalent  resistance 
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Sp  particle  surface  area  per  unit  volume 

T  effective  fuel  particle  temperature 
Tg  bulk  temperature  of  coolant 
Ts  fuel  particle  surface  temperature 
U  effective  over  all  heat  transfer  coefficient 
|i  viscosity 
v0  superficial  velocity 

V cv  size  of  the  control  volume 
W  mass  flowrate 

ACRONYMS 

ALMCR  Advanced  Liquid  Metal  Cooled  Reactor 
BNL  Brookhaven  National  Laboratory 
HTGR  High  Temperature  Gas  Cooled  Reactor 
IMEO  Initial  Mass  in  Earth  Orbit 
INEL  Idaho  National  Engineering  Laboratory 
LEO  Low  Earth  Orbit 

NERVA  Nuclear  Engine  for  Rocket  Vehicle  Application 
NTR  Nuclear  Thermal  Rocket 
PBR  Particle  Bed  Reactor 
PFNTR  Pressure  Fed  Nuclear  Thermal  Rocket 
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SNAP  Space  Nuclear  Applications  Program 
SNL  Sandia  National  Laboratory 
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